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Summary 


This  project  sets  out  to  identify  and  address  the  need  for  adaptive  ascent  guidance  tech¬ 
niques  necessary  for  responsive  launch.  The  required  short  response  time  and  potential  un¬ 
availability  of  mission  scenario  specifics  dictate  that  the  conventional  launch  ascent  guidance 
technology,  which  relies  heavily  on  pre-mission  planning  and  requires  significant  lead  time, 
is  simply  incompatible  to  the  needs  in  responsive  launch.  This  report  provides  comprehen¬ 
sive  details  to  two  recent  advanced  ascent  guidance  algorithms,  tailored  to  endo-atmospheric 
and  exo-atmospheric  optimal  ascent,  respectively.  These  algorithms  generate  optimal  ascent 
guidance  commands  based  on  the  current  state  of  the  vehicle,  the  currently  selected  target¬ 
ing  condition,  and  available  vehicle/environment  modeling  and  wind  information.  Should 
any  in-flight  changes  occur  in  vehicle  condition,  vehicle  health  status,  and  mission  objective 
(such  as  call  back  and  change  of  targeting  condition),  the  guidance  algorithms  would  be 
able  to  adapt  to  the  changes.  Throughout  the  report,  verification,  validation  and  extensive 
testing  of  the  algorithms  have  been  performed  with  many  mission  scenarios  and  a  number  of 
different  launch  vehicle  configurations,  including  winged  reusable  and  conventional  expand¬ 
able,  single-stage  and  multiple-stage,  launch  vehicles.  This  work  has  clearly  demonstrated 
that  promising  techniques  and  algorithms  have  reached  a  level  where  they  could  potentially 
lead  to  fully  automated  closed-loop  optimal  ascent  guidance  from  liftoff  to  orbital  insertion. 
It  is  strongly  recommended  that  continuing  efforts  be  made  to  further  enhance,  develop, 
and  mature  these  algorithms  and  techniques  to  pave  the  way  to  their  eventual  adaptation  in 
launch  operations. 
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Introduction 

2.1  State  of  Current  Launch  Ascent  Guidance  Tech¬ 
nology 

The  Air  Force’s  interests  in  achieving  operationally  responsive  launch  and  executing  time- 
critical,  global-reach  missions  from  space  require  far  greater  autonomy,  flexibility,  and  capa¬ 
bility  of  the  guidance  systems  than  currently  exist  for  launch  vehicles  and  entry  spacecraft. 
The  driving  motivation  for  this  research  is  that  the  challenges  for  realizing  responsive  access 
to  and  from  space  lie  not  only  in  hardware  and  operations,  but  also  equally  in  software 
and  algorithms.  Traditionally,  launch  and  entry  guidance  and  control  (G&C)  software  and 
parameters  are  designed  for  a  specified  mission,  payload,  and  targeting  condition.  In  ascent 
guidance,  the  current  technology  employs  open-loop  ascent  guidance  during  the  atmospheric 
portion  of  the  flight.  The  guidance  commands  are  generated  on  the  ground  based  on  the  or¬ 
bital  insertion  conditions,  vehicle  Modeling,  and  vehicle  load  and  integrity  constraints.  Prior 
to  launch,  the  day-of-launch  wind  profile  is  used  to  update  the  ascent  guidance  commands 
which  are  then  uploaded  into  the  launch  vehicle.  Such  elaborate  planning  and  updates  are 
essential  for  the  open-loop  guidance  to  ensure  that  the  load  limits  will  not  be  exceeded  dur¬ 
ing  ascent  in  the  inevitable  presence  of  winds,  the  performance  is  optimal,  and  the  targeting 
conditions  can  be  met.  But  this  is  a  time-consuming  and  labor-intensive  process,  done  well  in 
advance  of  the  actual  mission.  In  2004  the  Missile  Defense  Agency  reported  that  it  typically 
take  up  to  6  months  to  update,  re-validate,  and  check  out  the  guidance  and  control  software 
and  I-loads  when  targeting  conditions  and  mission  parameters  change.  See  MDA04-103, 
Research  Project  Call  “Flexible,  Rapid  Launch  Vehicle  Control  Software  Generation  and 
Checkout”.  The  following  is  the  excerpt: 

Current  inability  to  rapidly  generate  and  test  launch  vehicle  control  software  limits  MDA ’s 
capability  for  responding  to  late- developing  modifications  to  test  requirements  that  affect  bal¬ 
listic  missile  target  presentations.  If  modifications  to  test  requirements  involve  changes  to 
target  flight  dynamics,  trajectory,  separation  events,  etc.,  software  for  launch  vehicle  con¬ 
trol  must  be  re-generated,  verified,  re-installed,  and  validated  in  the  launch  vehicle  through 
check-out  procedures  to  include  integrated  testing.  For  multi-stage  vehicles,  typical  times  for 
developing  guidance,  navigation  and  control  (GN&C)  software  alone  approach  six  months... 
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In  crewed  launch  missions,  the  lead  time  required  for  planning  is  even  longer,  up  to  two 
years.  Below  is  the  personal  communication  from  a  senior  engineer  working  in  the  Flight 
Dynamics  Department  for  the  Space  Shuttle,  in  response  to  the  inquiry  on  Shuttle  launch 
planning  related  to  GN&C  and  the  man-power  requirement: 

The  Space  Shuttle  flight  design  process  can  take  anywhere  from  9  months  to  2  years, 
depending  on  the  mission.  Flights  that  involve  unique  payloads  are  more  complex,  the  most 
recent  example  that  I  can  think  of  would  be  STS-93  the  deploy  of  the  Chandra  observatory. 
That  payload  was  very  heavy  with  an  aft  CG  that  made  the  design  of  Ascent  Abort  trajec¬ 
tories  a  challenge,  particularly  the  glided  side  of  the  abort  with  full  OMS  pods,  and  a  heavy 
payload.  Recently  however,  all  the  flights  are  going  to  the  International  Space  Station  which 
has  simplified  the  flight  design  process  since  all  the  flights  now  have  similar  weight  and  tra¬ 
jectory  constraints  (51.6  deg.  Inclination,  ground  up  rendezvous  with  ISS,  and  medium  sized 
payloads  due  to  limitations  on  payload  mass  to  high  inclination  orbits).  This  has  also  helped 
with  the  flight  design  for  STS-300,  which  is  a  Launch  on  Need  (LON)  rescue  flight  to  the 
ISS.  If  the  orbiter  gets  damaged  and  can  not  return,  we  would  get  out  the  STS-300  trajectory 
products  and  refine  them  for  the  current  situation  and  launch  within  40  days. 

On  the  days  prior  to  launch  a  small  group  of  people  work  to  update,  many  of  the  ascent  I- 
loads  (see  below  for  a  definition  of  I-Load)  for  the  conditions  on  launch  day,  this  is  called  the 
” Day  of  Launch  I-Load  Update  (DOLILU,  pronounced  ’’doll  ee  lou”).  Atmospheric  conditions 
(winds,  and  atmospheric  density),  and  final  estimates  of  propellant  loading  all  factor  into 
these  final  adjustments  to  the  Hoads.  Items  that  get  updated:  First  stage  steering  commands, 
Throttle  bucket,  OMS  Assist  quantity,  etc.  DOLILU  team  6-10  people. 

The  Flight  Design  and  Dynamics  department  has  about  350  employees  that  do  everything 
from  simulation  software  modification  and  maintenance  to  the  actual  I-load  design.  Many  of 
these  people  are  dedicated  to  areas  that  would  not  apply  to  other  launch  vehicles  like  Ascent 
Aborts,  Entry,  and  Orbit  Operations.  Additionally,  there  is  a  significant  amount  of  effort 
that  goes  into  producing  crew  procedures,  and  many  of  the  people  in  this  department  also 
work  on  console  during  the  mission.  Also  they  are  staffed  to  support  a  flight  rate  of  6-8 
flights  a  year  plus  development  of  new  programs. 


2.2  Benefits  and  Objectives  of  the  Research 

Until  the  technical  challenges  in  update  and  design  of  G&C  algorithms  and  software  on 
a  short  notice  are  satisfactorily  addressed,  on-demand  launch  would  not  be  realistic  even 
for  a  vehicle  already  on  the  launch  pad  or  in  orbit.  Adaptive  closed-loop  ascent  guidance 
technology  through  the  atmosphere  is  not  only  desirable,  but  may  be  necessary  to  truly 
achieve  on-demand  responsive  space  launch. 

Listed  below  are  some  of  the  major  benefits  of  such  an  adaptive  closed-loop  ascent  guid¬ 
ance  technology  that  are  unmatched  by  the  open-loop  guidance,  and  directly  tied  to  the  Air 
Force’s  goal  of  responsive,  reliable  and  affordable  access  to  space: 

•  Significantly  shortened  pre-launch  guidance  preparation.  With  closed-loop  ascent  guid- 
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ance,  there  is  little  or  no  need  for  off-line  guidance  planning  and  analysis.  The  day- 
of-launch  wind  data  can  be  incorporated  into  the  guidance  algorithm.  Therefore,  pre¬ 
launch  guidance  command  update  is  unnecessary.  In  fact,  when  and  if  in-flight  means 
for  real-time  wind  speed  measurement  becomes  available,  the  measured  wind  profile 
could  be  directly  fed  to  the  closed-loop  ascent  guidance  algorithm  which  in  turn  would 
generate  corresponding  optimal  guidance  commands.  In  such  a  scenario  the  pre-launch 
process  of  taking  wind  measurements,  which  typically  takes  hours,  could  even  be  elim¬ 
inated. 

•  Operational  flexibility.  Closed-loop  ascent  guidance  would  be  fully  automated  and 
require  no  labor-intensive  pre-mission  analysis  and  re-planning  whenever  the  mission 
profile  changes.  Last-minute  changes  of  the  target  orbit  in  time-critical  missions  could 
be  easily  handled  by  the  guidance  system.  Even  in-flight  change  of  target  orbit  could 
be  possible. 

•  Dramatically  reduced  reoccurring  costs  related  to  guidance.  The  same  features  of  the 
closed-loop  ascent  guidance  that  provide  operational  flexibility  also  result  in  greatly 
reduced  need  for  human  intervention.  Thus  the  operational  costs  related  to  ascent 
guidance  could  be  reduced  to  minimum. 

•  Fault  Tolerance.  A  closed-loop  ascent  guidance  system  is  capable  of  adapting  to  se¬ 
vere  off-nominal  conditions.  It  could  readily  make  use  of  vehicle  health  information, 
accommodate  recoverable  system  failures,  and  still  successfully  complete  the  mission. 
An  example  is  the  failure  of  one  engine  (or  multiple  engines,  as  long  as  the  remaining 
engines  can  still  put  the  vehicle  into  orbit).  The  closed- loop  ascent  guidance  algorithm 
could  compute  new  optimal  guidance  commands  based  on  the  remaining  thrust  of  the 
vehicle,  and  guide  the  vehicle  to  fly  a  different  (longer)  trajectory  to  the  target  orbit. 

In  pursuit  of  the  required  technical  advances  in  realizing  adaptive  ascent  guidance,  the 
chief  goals  of  this  project  are  as  follows: 

1.  Identify  and  evaluate  advanced  algorithms  that  may  be  the  candidate  for  rapid  space 
launch  mission  planning  and  potentially  closed-loop  endo-atmospheric  ascent  guidance. 
Develop  new  algorithm  or  enhancements  to  an  existing  algorithm  where  appropriate. 

2.  Develop  robust,  reliable  and  fast  exo-atmospheric  ascent  guidance  algorithm  for  opti¬ 
mal  ascent  of  multi-stage  launch  vehicles.  The  algorithm  should  allow  full  autonomous 
guidance  for  optimal  coasts  between  two  burns/stages. 

3.  Integrate  the  endo-  and  exo-atmospheric  optimal  ascent  guidance  algorithms  for  com¬ 
plete  end-to-end  launch  operations. 

4.  Develop  testing  scenarios  and  evaluation  metrics  to  test  the  ascent  guidance  algorithms; 
validate  the  algorithms  with  other  existing  standard  aerospace  industry  trajectory 
software. 
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5.  Demonstrate  the  applicability  and  versatility  of  the  ascent  guidance  algorithms  to 
widely  different  launch  vehicle  configurations,  including  winged  reusable  launch  vehicle 
and  axial-symmetric  expendable  launch  vehicle. 

This  final  report  is  organized  in  accordance  with  the  above  stated  goals. 
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Coordinate  Systems 


The  need  to  define  appropriate  coordinate  systems  arises  from  two  considerations.  First, 
there  may  be  some  particular  coordinate  system  in  which  the  position  and  velocity  of  the 
flight  vehicle  “make  sense”.  For  navigation  we  are  concerned  with  position  and  velocity 
with  respect  to  the  Earth,  whereas  for  vehicle  aerodynamic  performance  we  need  position 
and  velocity  with  respect  to  the  atmosphere.  Second,  there  are  coordinate  systems  in  which 
the  phenomena  of  interest  are  most  naturally  expressed.  The  direction  of  a  rocket’s  thrust 
may  often  be  considered  fixed  with  respect  to  the  body  of  the  vehicle.  There  are  three 

coordinate  systems  used  in  this  report  -  Earth  Centered  Inertial  coordinate  system  (ECI), 
guidance  coordinate  system  (or  launch  plumbline  frame),  and  vehicle  body  frame.  This 
section  provides  descriptions  of  above  three  coordinate  systems  and  of  the  algorithms  used 
to  transform  quantities  between  different  systems. 


3.1  Earth  Centered  Inertial  Coordinate  System  XjYjZj 

As  its  name  suggests  this  coordinate  system  has  its  origin  at  the  center  of  the  Earth.  The 
z-axis  Zj  is  parallel  to  the  Earth’s  rotation  axis  (positive  to  the  North).  It  is  assumed  that 
the  x-axis  Xj  points  toward  the  intersection  of  the  Equator  and  the  Greenwich  Meridian 
at  the  time  of  launch  (this  is  a  simplification  that  can  be  easily  removed).  The  y-axis  Yj 
completes  the  right-hand  coordinate  system.  Thus  it  is  convenient  for  specifying  the  location 
of  ground  stations  and  ground-based  experiments  as  these  are  fixed  quantities  in  the  ECI 
system. 

When  ECI  coordinates  are  expressed  in  spherical  form,  the  latitude  component  is  identical 
to  what  is  termed  geocentric  latitude  by  astronomers  and  geographers.  However,  note  that 
this  is  different  to  the  system  of  geodetic  latitude  used  in  normal  map-making.  The  geodetic 
latitude  at  any  location  is  the  angle  between  the  equatorial  plane  and  the  local  normal  to 
the  Earth’s  surface.  In  general  that  normal  is  not  parallel  to  a  radius  vector  because  the 
shape  of  the  Earth  is  an  oblate  spheroid  and  not  a  sphere. 
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meridional  plane 


Figure  1:  Earth  centered  inertial  and  launch  plumbline  coordinate  systems 

3.2  Launch  Guidance  (Plumbline)  System  XpYpZp 


The  guidance  system  is  an  inertial  system,  which  is  also  called  the  inertial  launch  plumbline 
coordinate  system  whose  origin  is  at  the  center  of  the  Earth.  The  x-axis  XP  is  defined  from 
the  center  of  the  Earth,  parallel  to  the  gravity  direction  at  the  launch  site  and  positive  up 
(the  same  as  Xq  in  Fig.  1).  The  z-axis  Zp  is  pointing  downrange  along  the  launch  azimuth 
direction  and  the  y-axis  Yp  completes  the  right-hand  system.  The  XpYpZp  frame  has  its 
axes  parallel  (or  coincide  in  the  case  of  the  x-axis)  with  those  of  the  XqYqZq  frame  in  Fig.  1, 
except  that  the  latter  has  its  origin  at  the  launch  site.  The  longitude  and  geocentric  latitude 
of  the  launch  site  is  dehned  by  (@,<hc).  The  launch  azimuth  Az  for  an  ascending  orbit  is 
defined  by 


Az  =  sin 


-l 


cos? 

cos 


A,  = 


TT 


sm 


-l 


cos? 

cos 


for  ascending  orbit 

for  descending  orbit 


(1) 


where  ?  is  the  target  orbital  inclination. 
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North 


Figure  2:  Geodetic  and  geocentric  latitude 

Note  that  geocentric  latitude  <FC  is  used  instead  of  geodetic  latitude  The  geocentric 
latitude  of  a  point  is  the  angle  between  the  equatorial  plane  and  a  ray  through  the  point 
from  the  Earth’s  center.  The  geodetic  latitude  is  the  angle  between  the  local  zenith  and  the 
equatorial  plane.  Due  to  the  Earth’s  oblateness,  geodetic  latitudes  (the  most  common  form 
of  Earth  location)  are  slightly  greater  than  geocentric  latitudes  except  at  the  equator  and 
poles  where  they  are  identical.  The  relationship  between  <f>c  and  is  given  by 

tan  <f>c  =  (1  —  e2)  tan  (2) 

where  e  =  0.0818191  is  the  eccentricity  of  the  Earth. 


3.3  North-East-Down  (NED)  System 

This  is  an  inertial  system  with  the  origin  at  the  launch  site  at  the  launch  time.  The  x-axis 
points  to  the  North,  the  y-axis  to  the  local  East,  and  the  z-axis  completes  a  right-hand 
system  (down). 
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Figure  3:  Vehicle  body  coordinate  system  showing  Euler  angles 

3.4  Vehicle  Body  Coordinate  System  XbYbZb 

The  vehicle  body  coordinate  system  is  fixed  to  the  vehicle  as  shown  in  Fig.  3.  The  x-axis  Xp 
coincides  with  the  vehicle  body  longitudinal  axis.  The  z-axis  Zb  lies  in  the  plane  of  symmetry 
(or  some  reference  plane  in  the  case  of  asymmetric  shapes),  pointing  “downward” .  The  y-axis 
Yb  is  perpendicular  to  these  axes  forming  a  right-handed  coordinate  system.  Positive  Yb, 
thus  points  to  the  right  when  looking  forward.  The  origin  is  generally  taken  at  the  vehicle 
center  of  gravity  or  at  a  fixed  reference  location  relative  to  the  geometry. 

The  Euler  angles  are  also  shown  in  Fig.  3.  The  yaw,  pitch  and  roll  angles  are  denoted  by 
9,  -0,  and  0,  respectively. 


3.5  Coordinate  Transformation 

Let  Tep  be  the  coordinate  transformation  matrix  from  ECI  coordinate  system  to  the  plumbline 
frame. 


Tep  — 


cos  0  cos  $  sin  0  cos  $  sin  $ 

—  sin  0  cos  Az  +  cos  0  sin  $  sin  Az  cos  0  cos  Az  +  sin  0  sin  $  sin  Az  —  cos  $  sin  Az 

—  sin  0  sin  Az  —  cos  0  sin  $  cos  Az  cos  0  sin  Az  —  sin  0  sin  $  cos  Az  cos  $  cos  Az 
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(3) 


where  0  is  the  launch  site  longitude,  <h  is  the  launch  site  geocentric  latitude,  and  Az  is  the 
launch  azimuth  (see  Equation  1).  Using  Tpp,  we  can  easily  do  the  coordinate  transformation 
from  ECI  system  to  plumbline  frame,  and  vice  versa.  Using  rotation  sequence  of  pitch- yaw- 

roll  (also  referred  to  as  2-3-1  rotation),  the  coordinate  transformation  matrix  Tpp,  from  body 
frame  to  plumbline  frame  is 

cos  9  cos  ip  sin  9  sin  <p  —  cos  9  sin  ip  cos  <p  sin  9  cos  <p  +  cos  9  sin  ip  sin  <p 
Tp p  =  sin  ip  cos  ip  cos  <p  —  cos  ip  sin  <p  (4) 

—  sin  9  cos  ip  cos  9  sin  <p  +  sin  9  sin  ip  cos  <p  cos  9  cos  cp  —  sin  9  sin  ip  sin  cp 

where  9 ,  ip,  <p  are  the  three  Euler  angles  -  pitch,  yaw,  and  roll,  respectively.  The  three 
columns  of  Tpp,  are  indeed  the  three  unit  vectors  of  the  body  axes  1  ly,  lz  in  plumbline 
frame.  Therefore,  the  unit  vector  of  the  body  x-axis  is  given  by 

cos  6  cos  ip 

lb  =  sin-0  (5) 

—  sin  9  cos  ip 

The  unit  vector  of  the  body  y-axis  in  plumbline  frame  is  defined  as 

sin  9  sin  <p  —  cos  9  sin  ip  cos  <p 

ly  =  cos  ip  cos  cp  (6) 

cos  9  sin  cp  +  sin  9  sin  ip  cos  <p 

The  body  z-axis  can  be  determined  by  imposing  the  right-hand  rule 

lz  =  lb  X  ly  =  -In  (7) 

Once  we  find  the  three  body  axes  in  plumbline  frame,  the  Euler  angles  could  be  easily 
calculated  using  the  following  relationship. 


(8) 


where,  1  &x,  1  by  and  1  are  the  three  components  of  the  unit  vector  1&.  And  lyy,  lzy ,  are  the 
y-components  of  unit  vector  ly  and  lz  respectively. 
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Endo-Atmospheric  Closed-Loop 
Ascent  Guidance 

4.1  Introduction 

The  ascent  guidance  system  of  a  rocket-powered  launch  vehicle  determines  the  attitude 
commands  and,  when  applicable,  engine  throttle  command  during  the  ascent  of  the  vehicle. 
It  is  well  known  that  whether  or  not  the  ascent  trajectory  is  optimal  can  have  a  significant 
impact  on  propellant  usage  for  a  given  payload,  or  on  payload  weight  for  the  same  gross 
vehicle  weight.  Consequently  ascent  guidance  commands  are  usually  optimized  in  some 
fashion.  In  fact,  ascent  guidance  is  one  of  the  most  notable  engineering  fields  where  optimal 
control  theory  has  found  routine  applications.  Successful  vacuum  rocket  guidance  software 
based  on  the  optimal  control  theory  includes  the  Iterative  Guidance  Mode  (IGM)  for  the 
Saturn  rockets,1  and  the  Powered  Explicit  Guidance  (PEG)  for  the  Space  Shuttle.2  These 
algorithms  solve  the  optimal  vacuum  powered  flight  problem  on-board  in  each  guidance 
update  cycle  using  the  current  condition  as  the  initial  condition  of  the  solution.  Therefore 
the  guidance  strategy  in  effect  is  closed-loop. 

One  of  the  major  open  challenges  of  ascent  guidance  lies  in  the  endo-atmospheric  por¬ 
tion  of  the  flight.  The  presence  of  the  aerodynamic  forces,  loads  and  winds  significantly 
complicates  the  optimal  ascent  problem,  making  the  solution  process  much  more  difficult  to 
converge  reliably  and  sufficiently  fast  for  on-board  applications.  For  these  reasons  typical 
current  ascent  guidance  inside  the  atmosphere  is  open-loop.3  In  such  an  approach,  the  guid¬ 
ance  commands  are  generated  off-line,  updated  with  the  day-of-launch  wind  data  prior  to 
launch,  and  loaded  into  the  launch  vehicle  for  use  during  the  ascent  through  the  atmosphere. 
While  very  successful  in  nominal  ascent  guidance,  the  open-loop  approach  inherently  lacks 
the  adaptive  capability  to  handle  contingencies  and  aborts,  even  with  extensive  off-line  plan¬ 
ning  at  great  costs.  Open-loop  guidance  also  does  not  possess  the  robustness  necessary  to 
cope  with  significant  off-nominal  conditions  and  system  Modeling  uncertainty,  especially  for 
new  launch  vehicles  for  which  little  or  no  flight  data  is  available.  The  required  re-planning 
and  re-generation  of  the  open-loop  ascent  guidance  commands  whenever  any  mission  or 
system  parameters  change  are  costly  in  both  developmental  and  operational  phases  of  the 
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launch  vehicle. 

A  closed-loop  ascent  guidance  algorithm  could  address  all  of  the  above  deficiencies  of 
open-loop  guidance.  The  search  for  a  feasible  algorithm  to  solve  the  optimal  control  prob¬ 
lem  on-board  for  closed-loop  atmospheric  ascent  guidance  dates  back  to  the  1960s.  The  work 
by  Brown  and  Johnson  represents  one  of  the  earliest  attempts  in  this  direction.4  In  a  series 
of  recent  work5-'  that  has  stimulated  renewed  interest  in  this  area,  Calise  et  al  develop  a 
hybrid  approach  to  the  problem.  In  this  approach  the  analytical  solution  of  the  optimal 
vacuum  flight  and  numerical  collocation  for  atmospheric  portion  are  combined.  The  vacuum 
solution  serves  as  the  initial  guess  for  the  atmospheric  flight,  and  a  homotopy  method  is 
used  to  gradually  phase  in  the  aerodynamic  terms  and  path  constraint-related  terms.  Ref¬ 
erences8-11  contain  further  recent  enhancements  and  more  development  in  endo-atmospheric 
ascent  guidance. 

This  part  of  the  report  describes  the  endo-atmospheric  guidance  algorithm  used  in  this 
work.12  The  main  components  of  this  part  of  the  report  are: 

•  Comprehensive  treatment  to  the  3-dimensional  optimal  ascent  problem  subject  to 
the  common  path  constraints  and  orbital  insertion  conditions,  and  techniques  to  ad¬ 
dress  a  number  of  on-board  implementation  issues,  some  of  which  are  unique  to  non- 
axisymmetric  launch  vehicles  (such  as  a  lifting-entry  reusable  launch  vehicle); 

•  Demonstration  of  the  suitability  for  the  solution  of  the  ascent  guidance  problem  by 
a  classical  finite-difference  approach  which  can  be  interpreted  as  a  special  form  of 
collocation,  but  is  conceptually  simpler  and  easier  to  implement; 

•  Illustration  of  the  capability  and  feasibility  of  on-board  closed-loop  ascent  guidance  by 
a  series  of  carefully  designed  tests. 


4.2  Ascent  Trajectory  Dynamics 


The  equations  of  motion  for  a  rocket-powered  launch  vehicle,  in  a  central  gravitational  field, 
expressed  in  an  inertial  coordinate  system  are  as  follows: 


r 

=  V 

(1) 

V 

,  v  Tlb  A  N 

=  q(r)  +  +  + 

m(t)  m(t )  m(t ) 

(2) 

tin 

9^ vac 

9oRp 

(3) 

where  r  and  V  E  R3  are  the  inertial  position  and  velocity  vectors;  g  the  gravitational 
acceleration;  Tvac  the  full  vacuum  thrust  magnitude;  rj  >  0  is  the  engine  throttle;  T  the 
current  thrust  magnitude  including  effects  of  throttle  modulation  and  thrust  loss  clue  to 
back  pressure.  In  this  formulation  the  total  engine  thrust  is  assumed  to  be  aligned  with  the 
body  longitudinal  axis,  and  is  not  gimbaled  independently.  The  vectors  A  and  N  are  the 
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aerodynamic  forces  in  the  body  longitudinal  and  normal  direction,  respectively;  1&  the  unit 
vector  defining  the  launch  vehicle  body  longitudinal  axis;  m(t)  is  the  mass  of  the  launch 
vehicle  at  the  current  time  t.  The  specific  impulse  of  the  engine  is  Isp  and  go  represents  the 
gravitational  acceleration  magnitude  on  the  surface  of  the  Earth. 


4.3  Definition  of  Vehicle  Body  Axis  Frame 


The  definition  of  vehicle  body  axis  frame  Xi/y^Zb  depends  on  how  we  want  the  vehicle  to  fly. 
We  can  choose  to  construct  lz  (therefore  the  symmetric  plane)  so  that  the  vehicle  flies  a 
zero  degree  (“heads-up”)  or  180  degree  ( “heads-down” )  bank  angle  trajectory. 


1 


y  ~ 


1  b~xr 

lb  x  r* || 


(4) 


The  Shuttle  adopts  this  heads-down  option.  This  heads-down  position  assists  in  communica¬ 
tions  with  the  ground  and  allows  instruments  within  the  cargo  bay  to  be  pointed  back  toward 
the  Earth,  which  is  required  for  many  of  the  experiments  carried  within  the  bay.  There  is 
probably  also  some  psychological  benefit  to  the  crew  since  they  are  given  spectacular  views 
of  home  rather  than  staring  into  the  cold  darkness  of  the  great  void  of  space. 

We  can  also  construct  lz  so  that  the  vehicle  flies  at  zero  sideslip  angle.13  The  definition 
of  vehicle  body  axis  frame  in  this  study  follows  this  zero-sideslip  formulation.  In  this  formu¬ 
lation,  the  launch  vehicle  symmetric  plane  is  assumed  to  be  always  the  plane  formed  by  the 
body-axis  1&  and  the  Earth  relative  velocity  vector  V  r.  Thus  the  sideslip  angle  remains  zero. 
Note  that  such  a  body-frame  orientation  necessitates  a  roll  angle  about  the  longitudinal  axis 
lb  to  null  the  sideslip  in  the  presence  of  cross  winds.  Physically,  this  is  the  so-called  “fly  into 
the  wind”  maneuver  as  in  the  case  of  Space  Shuttle  ascent. 


y- 


Thus  the  unit  vector  of  the  body  x-axis  is  the  same  as  1&;  the  unit  vector  of  the  body 
axis  is  defined  as 


1 


y  ~ 


lyr  x  1;, 
II  Ivy  x  lb 


(5) 


where,  lyr  =  Vr/Vr  is  the  unit  vector  in  the  direction  of  Vr.  The  unit  vector  of  the  body 
z-axis  completes  the  right-hand  system  lz  =  If,  x  ly.  Denote  the  body-normal  unit  vector 
by  ln  =  -lz.  Then 

(If,  X  Vr)  , 

1“  =  U  x  1  V  v  \\  (  a>0  )  (6) 

where,  V r  =  V  —  uJe  x  r  =  Earth  relative  velocity;  uJe  =  Earth  rotation  vector  in  the 
plumbline  launch  frame.  Note  in  this  formulation,  the  sideslip  (3  =  0.  Figure  4  shows  this 
configuration.  Clearly, 


cos  a  =  li  lVr  or  |sin  or|  =  ||ly_  x  lt 


(7) 


To  avoid  an  instantaneous  180-degree  rotation  of  ln, 
defined  to  be 
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n 


lb  X 


{Vr  X  lfe) 
||  lfo  X  Vr\\ 


when  a  crosses  a 
(  a  <  0  ) 


0,  ln  should  be 


(8) 
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Figure  4:  Launch  vehicle  body  frame  with  relative  velocity 


The  following  expression  for  ly  is  preferred  to  Equation  (5) 

ly  —  lVr  x  lft/ since  (9) 

The  reason  is  that  this  definition  is  valid  for  both  a  >  0  and  a  <  0  without  causing  the 
instantaneous  180-degree  rotation  in  ly,  when  1  yr  and  1&  cross  over  each  other  (a  changes 
sign) . 

4.4  Nondimensionalization 

For  better  numerical  conditioning,  the  following  nondimensionalization  is  used: 

•  The  distances  are  normalized  by  Rq,  the  radius  of  the  Earth  at  equator; 

•  Time  is  normalize  by  \J Ro/go ; 

•  The  velocities  are  normalized  by  y/Rogo,  the  circular  velocity  around  the  Earth  at  Ro 

The  gravity  is  Modeled  by  the  Newtonian  central  gravity  field.  With  some  abuse  of  notation, 
we  use  the  same  names  hereafter  for  the  dimensionless  variables.  Linder  the  zero-sideslip 
formulation  as  shown  in  Fig.  4,  the  dimensionless  equations  of  motion  from  Eqs  (1)  and  (2) 
become 

(  r'  =  V 

{  1  (10) 
\V'  =  -^  r  +  (T-A)lb  +  Nln 
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where  the  differentiation  is  with  respect  to  the  dimensionless  time.  Now  A  and  N  are  the 


aerodynamic  accelerations  in  g0  in  the  body  longitudinal  and  normal  direction  respectively, 
and  T  the  magnitude  of  the  thrust  acceleration  in  r/0 •  The  magnitudes  of  the  dimensionless 
aerodynamic  and  thrust  accelerations  are  given  by 

A  =  ^TPWV2S«/CJ(M0,a) 

(ii) 

N  =  ^r- p(r)Vr2SrefCN(Ma,a ) 

(12) 

T  =  [rjTvac  +  AT(r)\/m(t)g0 

(13) 

where,  p{r)  is  the  dimensional  atmospheric  density  at  radius  r; 
dimensionless  Earth-relative  velocity 

Vr  is  the  magnitude  of  the 

Vr  =  V  —  uje  x  r  —  Vw 

(14) 

where  Vw  is  the  wind  velocity  vector  and  ufo  is  the  Earth  angular  rotation  rate  vector.  The 
axial  and  normal  aerodynamic  coefficients  Ca  and  Cn  are  functions  of  Mach  number  Ma  and 
angle  of  attack  a.  They  are  expressed  in  analytical  forms  by  curving-htting  the  tabulated 
data.  The  thrust  loss  AT  due  to  the  back  pressure  is  a  function  of  altitude  through  the 
dependence  of  AT  on  ambient  pressure.  Note  that  the  mass  flow  rate  will  be  reduced  by  the 
same  percentage  as  the  thrust  when  the  thrust  is  throttled  down. 

4.5  Guidance  Problem  and  Constraints 

The  current  conditions  r0  and  Vo  are  assumed  to  be  known.  The  ascent  guidance  problem 
is  to  fold  the  desired  body-axis  orientation  l&(t)  at  each  instant  which  determines  the  thrust 
direction  and  aerodynamic  forces  during  the  atmospheric  portion  of  the  ascent.  The  engine 
throttle  rj  will  also  need  to  be  determined  for  enforcing  some  of  the  path  constraints,  as  will 
be  discussed  later.  The  final  conditions  will  be  the  engine-cutoff  conditions  which  ensure 
insertion  into  the  required  orbit.  These  orbital  insertion  conditions  can  in  general  be  written 
as  k,  0  <  k  <  6,  algebraic  end  conditions 

*(r(t/), V(t,))  =  0,  teS*  (15) 

Specifics  on  (15)  will  also  be  discussed  later.  In  addition,  there  will  be  path  constraints  on 
the  trajectory  for  safety  and  vehicle  integrity.  The  three  most  common  path  constraints  in 
ascent  guidance  will  be  considered:  the  product  of  dynamic  pressure  and  ck,  axial  thrust 
acceleration,  and  dynamic  pressure 

|  qoi\  <  Qa  (16) 

T  <  Tmax  (17) 

Q  —  Qmax  (IS) 

where  q  =  pVx  / 2.  The  constants  Qa,  Tmax ,  and  qmax  are  the  respective  limits  for  each  of  the 
corresponding  constraints.  Another  common  path  constraint  on  qfi | ,  the  product  of  dynamic 
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pressure  and  sideslip  angle,  is  not  included  here  because  in  our  zero-sideslip  formulation,  this 
quantity  usually  already  has  small  magnitude.  The  constraints  on  normal  acceleration  \N\ 
and  angle  of  attack  a,  if  necessary,  can  be  handled  in  the  same  way  as  constraint  (16)  is 
handled.  Thus  they  will  not  be  discussed  separately.  Collectively,  the  above  path  constraints 
may  be  written  in  a  compact  form 

S(r,  V,  lb,t)  <  0  (19) 

4.6  Optimal  Control  Problem  and  Necessary  Condi¬ 
tions 

The  mathematical  tool  used  to  find  the  optimal  ascent  guidance  commands  is  the  optimal 
control  theory.  In  this  setting  a  performance  index  is  defined.  The  minimization  of  this 
performance  index  is  usually  tied  in  one  way  or  the  other  to  the  minimization  of  propellant 
usage.  Denote  the  performance  index  by 

J  =  </>(rf,Vf,tf)  (20) 

where  tf  is  the  engine  cutoff  time  and  r f  and  V /  are  the  position  and  inertial  velocity  of 
the  launch  vehicle  at  tf.  The  functional  form  of  J  is  best  selected  to  be  most  convenient  for 
a  particular  formulation  of  the  optimal  ascent  problem.  Typical  choices  are  J  =  tf,  for  the 
minimum-time  problem,  or  J  =  1  /r f  —  Vj .  for  the  maximum-energy  problem  with  a  fixed 

tf- 

Most  recently  it  has  been  shown  that  in  general  no  singular  optimal  thrust  programs  exist 
in  atmospheric  ascent.8  Therefore  the  optimal  throttle  is  bang-bang  type.  In  this  report  the 
engine  throttle  r/  is  treated  as  a  given  (possibly  time- varying)  input.  Thus  the  variation  of 
the  mass  m{t)  is  considered  a  prescribed  function  of  time,  not  a  state.  In  order  to  focus 
on  the  presentation  of  the  essentials  of  the  approach,  the  path  constraints  Eqs.  (16-18)  will 
be  added  later.  With  these  assumptions  and  noting  the  constraint  1^1 b  =  1,  we  define  the 
Hamiltonian 

H  =  p^V  +  py  +  (T  —  A)lb  +  Nln  +  p( l^lfe  -  1)  (21) 

where  p  is  a  scalar  multiplier,  and  pr  and  pv  G  R3  are  the  so-called  costate  vectors.  Let 
the  asterisk  signify  the  optimal  values  of  the  relevant  variables.  The  standard  necessary 
conditions  for  the  optimal  solution  are14  (using  the  notion  of  Maximum  Principle) 

Pr 

Pv 

H(pr,Py,r*,V*,ll,t) 


dH 

dr 

dH 

~dV 

=  m&xH(pr,pv,r*,V*,lb,t) 

1  b 
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(22) 

(23) 

(24) 


The  derivation  of  the  expressions  of  the  costate  equations  (22-23)  is  quite  involved.  The 
detailed  equations  are  provided  in  Appendix  A.  The  optimal  solution  must  also  satisfy  the 
terminal  constraints  (15)  and  the  following  transversality  conditions 


Pr(tf) 

Pv(tf ) 


drf 

d<t>(r  f,  Vf,  tf) 


dV 


+ 


+ 


d<t>(rf,  V/,tf)  ,  ( <9T 
drf 
<9T 
3VL 


/ 


v 


H(pr,pv,r*,V*,rb,t)\tf  =  ^- 


d(j) 
f 


(25) 

(26) 
(27) 


where  v  G  Rk  is  a  constant  multiplier  vector.  The  last  condition  (27)  is  for  the  cases  where 
the  final  time  tf  is  not  specified.  The  first  two  conditions  (25-26)  can  be  combined  to 
eliminate  the  unknown  vector  v  and  yield  6  —  k  independent  conditions  involving  only  final 
costate  Pf  =  {pj.f  Pvf)T  and  final  state  Xf  =  (rj  Vj)T .  The  general  approach  will  be  first 
finding  the  6  —  k  linear  independent  solutions  of  the  homogeneous  system 


<9T 
d  xf 


£  =  0 


Let  £i(xf)  G  R6,  i  =  1,...,6  —  k  be  such  solutions.  Note  that  £j’s  are  functions  of  Xf. 
Transversality  conditions  (25)  and  (26)  are  then  equivalent  to 


(pf  +  &  =  ri(P/,  ^/)  =  0,  i  =  1,  •  •  • ,  6  —  k.  (28) 

The  k  terminal  constraints  (15)  plus  the  above  6  —  k  conditions  constitute  the  6  terminal 
conditions  for  the  optimal  control  problem.  For  a  given  problem,  the  conditions  in  (28)  can 
often  times  be  obtained  more  conveniently  by  using  the  terminal  constraints  (15)  and  taking 
dot  products  of  Eqs.  (25)  and  (26)  with  appropriate  vectors  related  to  the  final  state  Xf. 
Examples  will  be  given  later. 


The  optimality  condition  (24)  necessitates 


dH 
<91  b 


=  0 


(29) 


Define  s  —  ||  If,  x  lyr||.  The  expansion  of  the  optimality  condition  dH/dlb  =  0  requires, 
among  other  relationships,  the  following 


da 

1  b 

dln 

<91  b 


cos  a 


sm  a 


sm  a 


Lvr 


1 


~  1  (lyrlf>)l3  +  Ifcly,  +  —  [(lyrl&)l&  “  Vr]  [(lyrlfc)lvr 


(30) 

(31) 


where  I3  is  a  3  x  3  identity  matrix.  Let 


a  =  pv[(lvrlb)(lpvlb)  ~  (lvrlPy)]/s 
b  =  -pvAa  +  a,Na 
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where  lpv  =  pv/pv,  Aa  =  dA/da  and  Na  =  dN/da.  When  evaluating  dH/dlb,  keep  in 
mind  that  both  A  and  N  are  functions  of  a,  therefore  functions  of  1*,.  Carrying  out  the 
differentiations  and  collecting  terms,  we  eventually  have 


1  * 

1b  ~ 


{-[T-A  +  N(lrvVb)/s\pv 


(2 n  +  b/  tana  —  aN/s 2) 

+  [b/  sin  a  -  Npv{lTpvTb)/ s  -  aN (l^lj). / s2j  lVr} 

=  C1(x,p,rb)pV+C2(x,p,ll)Vr 


(32) 


where  c\  and  c2  are  scalar  functions  of  the  state,  costate,  and  l*b.  Hence  we  conclude  that 
the  optimal  body-axis  lies  in  the  plane  formed  by  the  primer  vector  pv  and  relative  velocity 
vector  Vr.  A  similar  conclusion  is  reached  by  using  a  geometric  approach  in  an  earlier 
work  by  Vinh,15  where  the  thrust  direction  and  aerodynamic  force  vector  are  assumed  to  be 
independent  controls. 

The  condition  (32)  suggests  that  the  search  for  the  optimal  body-axis  orientation  can 
be  reduced  to  a  one-dimensional  search  in  the  plane  of  pv  and  Vr.‘  Let  <f>  be  the  angle 
between  the  vectors  pv  and  Vr.  At  each  instant  with  given  state  and  costate,  <f>  is  known. 
Denote  by  lPr  and  lpv  the  unit  vectors  in  the  directions  of  pr  and  pv,  respectively.  Notice 
that  Iff lpv-  =  cos(<3>  —  a)  and  1^1  pv  =  sin($  —  a).  Using  these  two  relationships  in  the 
expression  of  H  in  Eq.  (21),  it  is  clear  that  maximizing  H  with  respect  to  If,  is  equivalent 
to  dH/da  =  0,  which  in  turn  results  in' 

tan(<3>  —  a)(T  —  A  +  Na )  —  (Aa  +  N)  =  0  (33) 


Since  A,  N,  Aa  and  Na  are  generally  functions  of  a,  the  above  equation  needs  to  be  solved 
numerically  for  a. 

The  cases  for  canted  or  gimbaled  thrust  vector  are  not  included  in  above  formulation. 
But  the  same  conclusion  can  be  reached  in  a  similar  fashion  as  in  above  analysis.  This  result 
is  an  extension  of  the  well-known  primer  vector  theory  on  optimal  rocket  flight  in  vacuum.31 

Once  a  is  found  from  Eq.  (33),  c±  and  c2  in  Eq.  (32)  can  be  solved  in  terms  of  a  and  <3? 
by  taking  the  dot  product  of  (32)  with  lpv  and  with  lVr 


1*  — 
1b  ~ 


/  sm  a 
\  sin  <f> 


lpv-  + 


sin  ( T  —  a) 
sin  <f> 


LVr 


(34) 


Note  that  as  the  atmospheric  density  decreases  (approaching  vacuum  flight),  the  aerody¬ 
namic  terms  diminish,  and  a  — >  <f>  from  Eq.  (33)  .  The  optimal  body  axis  in  (34)  and 
therefore  the  optimal  thrust  vector  become  aligned  with  the  primer  vector  pv.31 


4.7  Adding  the  Path  Constraints 

When  inequality  constraints  (19)  are  added  to  the  problem,  the  costate  equations  (22-23) 
will  have  additional  terms  related  to  the  constraints,  and  the  condition  to  determine  the 
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controls  will  change.  We  shall  discuss  each  of  the  three  constraints  in  Eqs.  (16-18)  in  the 
following. 

Constraint  S\  =  qa  —  Qa  <  0 


Without  loss  of  generality  the  absolute  sign  in  Eq.  (16)  is  removed  for  the  simplicity 
of  discussion.  This  constraint  is  a  zeroth  order  constraint  in  that  the  control  1;,  appears 
explicitly  (through  a)  in  the  constraint  itself.  In  such  a  case,  the  costate  equations  take  the 
form  of 


where  the  multiplier  Xqa  = 
Si  <  0.  When  Si  =  0,  Xqa 


P  = 


dH 


X  n 


dS  i 


(35) 


dx  qa  dx 

0  and  the  problem  is  the  same  as  in  the  preceding  sections  when 
satisfies  the  condition 


dH  |  x  dSx 

dlb  +  qadTb 


0 


(36) 


Using  condition  (30)  we  have 


dSi  da 

WT  =  =  q 

dlb  1  b 


cos  a 


sm  a 


(37) 


Using  condition  (36)  and  the  result  in  Eq.  (37)  and  following  the  steps  similar  to  those  in 
Section  3.6,  we  can  show  that  the  optimal  body  axis  lj)  is  again  in  the  plane  of  pv  and  Vr 
as  before.  Note  that  this  conclusion  cannot  be  assumed  to  be  always  true  for  any  inequality 
constraints,  and  must  be  verified  for  each  case.  In  this  case  condition  (36)  is  equivalent  to 


Or, 


dH  x  dS1  n 

- Xqa— —  —  0 

da  da 


Xqa 


dH/da 


where 


dH 

da 


=  pv[(T  —  A  +  Na)  sin($  —  a)  —  ( Aa  +  N)  cos($  —  a) 


In  a  finite  interval  where  S\  —  0,  the  angle  of  attack  is  directly  obtained  from 


(38) 


(39) 

(40) 


«  =  Qa/q  (41) 

The  direction  of  the  body  axis  If,  is  determined  by  a  in  the  plane  of  pv  and  Vr  by  Eq.  (34) 
as  before.  The  multiplier  Xqa  is  calculated  from  Eq.  (39)  and  used  in  the  costate  equation 
(16)  to  propagate  the  costate. 


Constraint  S2  =  T  —  Tmax  <  0 


Typically  for  the  usual  values  of  Tmax ,  this  constraint  only  becomes  active  after  the 
trajectory  is  outside  the  dense  atmosphere  where  A T(r)  ~  0.  Thus  whether  or  not  and 
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when  this  constraint  will  become  active  are  uniquely  determined  by  the  prescribed  engine 
throttle,  not  influenced  by  any  of  the  trajectory  state  variables  and  the  control  If,.  In  this 
sense  it  is  not  a  state  or  control  constraint  (recall  that  the  engine  throttle  is  regarded  as  a 
prescribed  input).  Once  S2  =  0,  the  engine  throttle  will  be  adjusted  according  to 

H  _  Tmaxrn(t)gQ 

T,!ac 

to  keep  S2— 0.  The  costate  equations  are  the  same  as  in  (22-23),  and  the  optimal  body  axis 
11  is  found  in  Eq.  (34)  as  before. 

Constraint  S3  =  q  —  qmax  <  0 


This  is  a  first-order  constraint  because  the  control  If,  appears  in  the  first-order  derivative 
of  S3.  Recall  that  q  =  p(r)V///2,  and  V r  is  defined  in  Eq.  (14).  So 


S'3  =  -prV?rTV  +  pVTrV'r 


(43) 


where  pr  =  dp/ dr  and 

V'r  =  V'-uEx  V-V'w  (44) 

Suppose  that  S3  —  0  in  a  finite  interval  [0,  t2\.  The  costate  equations  in  this  interval  are 


P  = 


dH  dS'3 
dx  q  dx 


The  corresponding  optimality  condition  is 


f  +  =  0 

dlh  dlb 


The  costate  will  have  a  jump  at  t\ 


P(K)  =p(t  1 )  +  K 


dS3 

dx 


where  n  is  a  constant  multiplier.  It  can  be  shown  that 


8S’ 

dh, 


d\(x ,  1  b)X^  r  T  d2(x ,  lf,)lf, 


(45) 


(46) 


(47) 


(48) 


where  d±  and  d2  are  two  scalar  functions.  Therefore  using  the  derivations  in  Section  3.6,  we 
can  show  that  condition  (46)  still  results  in  that  the  optimal  body  axis  lj(  lies  in  the  plane 
of  pv  and  Vr.  Conditions  similar  to  Eqs.  (38)  and  (39)  can  now  be  readily  derived.  And 
the  condition  S'3  =  0  provides  the  equation  for  the  solution  of  a  in  [fi,f2]- 

With  the  above  equations  laid  out,  the  optimal  solution  in  principle  can  be  found  nu¬ 
merically.  But  two  implementation  issues  arise:  1)  the  engine  throttle  is  known  to  be  more 
effective  in  regulating  the  dynamic  pressure  by  slowing  down  the  increase  of  the  velocity;  2) 
the  jump  condition  Eq.  (47)  makes  it  necessary  to  accurately  estimate  the  time  t\  in  the 
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solution  process  for  it  to  converge  quickly.  The  second  issue  may  not  be  a  problem  if  only 
nominal  ascent  is  considered,  because  a  good  estimate  of  t\  can  be  obtained  off-line  if  the 
vehicle  modeling,  day-of-launch  wind  data,  and  mission  parameters  are  well  known.  But  a 
chief  potential  advantage  of  closed-loop  ascent  guidance  is  for  autonomous  abort  guidance. 
In  aborts,  the  conditions  are  inevitably  well  off  nominal,  and  any  off-line  estimates  of  t\ 
could  be  no  better  than  an  arbitrary  guess. 

In  an  on-board  environment,  the  entire  ascent  trajectory  cc(-),  along  with  the  control 
lft(-)  and  throttle  rj(-),  is  generated  from  the  current  condition  x(t)  to  the  target  condition 
in  each  guidance  cycle.  The  current  attitude  commands  and  throttle  command  are  from 
1  b(t)  and  77(f),  the  first  data  point  in  the  guidance  solution.  To  address  the  above  issues  and 
keep  the  algorithms  reasonably  simple  and  robust  for  on-board  guidance  purposes,  we  adopt 
the  following  approach:  The  optimal  body-axes  in  the  guidance  solution  are  still  determined 
as  in  Sections  3.5  with  the  prescribed  throttle  where  no  constraint  on  the  dynamics  pressure 
is  considered.  This  is  justified  because  we  conclude  that  the  optimal  body  axis  remains  in 
the  plane  of  pv  and  Vr  even  when  q  <  qmax  is  active.  Next,  we  consider  the  first-order 
derivative  of  q  at  t  that  can  obtained  as 

q\t)  =  aq(x,  lb)q (t)  +  bq(x,  lb)  (49) 


where  77  is  the  current  engine  throttle,  and  the  expressions  for  aq  and  bq  are  readily  available 
when  one  makes  the  substitutions.  Let  <5  >  0  be  a  time- increment.  A  first-order  approxima¬ 
tion  of  q(t  +  5)  is 

q(t  +  ^)  «  q(t)  +  q'(t)5  =  q(t)  +  [aq(x,  lb)q (t)  +  bq(x,  lb)]<5  (50) 


When  deciding  the  throttle  command  at  the  current  time  t,  we  require  that  q(t  +  S)  <  qmax- 
Using  the  above  approximation  for  q{t  +  5)  in  this  condition,  we  have 


v{t)  < 


qmax  <?(^)  bq5 
aq5 


A 

=  Vq 


(51) 


Most  likely  77  will  have  a  minimum  allowable  setting  rjnnn  >  0  during  engine-on  period.  Let 
Vprb  >  Vmin  be  the  otherwise  prescribed  throttle  setting  (e.g.,  r)max )•  The  current  engine 
throttle  command  is  determined  by 


V 


Vprbi  If  Vq  V prb 

Vqi  if  Vmin  5;  Vq  —  V prb 
Vmin i  if  Vq  ^  Vmin 


(52) 


The  last  case  in  (52)  would  be  when  the  dynamic  pressure  constraint  cannot  be  met  by 
lowering  the  engine  throttle  within  the  allowable  range,  which  is  an  unlikely  event. 

Remarks: 


1.  Since  a  first-order  expansion  (50)  is  used  in  deriving  (51),  one  may  be  tempted  to 
suspect  that  the  enforcement  of  q  <  qmax  by  this  approach  is  acceptable  only  when 
5  is  small.  On  the  contrary,  it  can  be  shown  that  as  long  as  the  trajectory  begins  in 
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the  region  where  q  <  qmax ,  the  throttle  given  by  (52)  guarantees  strict  satisfaction  of 
the  constraint  q  <  qmax  for  any  5  >  0.  It  turns  out  that  when  q  — >  qmax,  V  =  Vq 
drives  q'  — >  0  for  any  5  >  0,  and  whenever  q  =  qmax ,  q'  =  0.  A  detailed  discussion  of 
strict  enforcement  of  first-order  constraints  by  this  type  of  technique  can  be  found  in 
a  previous  work  on  constrained  nonlinear  control  systems.16 

2.  Because  of  the  above  property,  an  appropriate  value  6  can  always  be  chosen  for  a  given 
vehicle  to  eliminate  the  undesirable  jitters  in  throttle  command  when  the  dynamic 
pressure  constraint  is  active,  and  still  accurately  enforce  the  constraint.  Typically  there 
is  a  minimum  value  of  6  above  which  a  6  will  render  the  engine  throttle  commands 
from  (52)  to  be  sufficiently  smooth. 

3.  The  benefits  of  this  approach  for  enforcing  the  constraint  q  <  qmax  are  its  simplicity, 
robustness,  and  effectiveness.  And  there  is  no  need  for  the  estimate  of  the  instant 
ti  when  the  constraint  becomes  active.  It  is  closed-loop  in  nature,  unlike  the  open- 
loop  “throttle  bucket”  approach  currently  in  use  for  the  Shuttle.  It  should  be  pointed 
out,  however,  that  the  adjustment  of  the  throttle  is  done  outside  the  solution  process 
of  the  optimal  control  problem.  Thus  the  guidance  solution  in  the  period  when  the 
dynamic  pressure  constraint  is  active  is  not  strictly  optimal  in  theory  (although  the 
body  axis  is  still  determined  according  to  the  necessary  conditions  of  the  optimal 
control  problem).  But  this  period  is  always  relatively  short  compared  to  the  total 
burn  time.  The  impact  of  this  short  deviation  from  the  theoretical  optimal  conditions 
on  performance  appears  to  be  negligible.  We  have  used  two  different  state-of-the-art 
trajectory  optimization  software  packages17, 18  to  cross-check  the  performance  of  the 
trajectories  obtained  herein.  The  differences  in  terms  of  propellant  consumption  have 
been  minimal  (See  numerical  results  later).  A  similar  observation  is  also  made  by 
Corvin19  on  the  effects  of  adjusting  throttle  to  control  dynamic  pressure. 


22 


5 


Numerical  Method 


5.1  Finite  Difference  Approach 

Finite  difference  is  one  of  several  classical  techniques  used  for  two-point-boundary- value 
problems  (TPBVPs).20  For  ascent  guidance  applications,  we  have  found  that  this  classical 
approach  is  well  suited  for  the  problem. 

After  substituting  the  control  If,  of  Eq.  (34)  in  the  equations  of  motion  and  costate 
equations,  and  denoting  y  =  ( xT  pT)T  G  R2n  with  n  =  6,  the  complete  TPBVP  is  now 

y'  =  f(t,y)  (i) 

B0(y0)  =  0  (2) 

Bf(yf)  =  o  (3) 

where  B o  =  x(to)  —  Xo  in  our  problem  represents  the  given  initial  condition  ,  and  B  f(y j)  =  0 
are  the  6  final  conditions  from  combining  the  orbital  insertion  conditions  Eq.  (15)  and  the 
transversality  conditions  (28).  Let  tf  be  the  specified  final  time.  The  TPBVP  is  to  find  a 
solution  y(t)  that  satisfies  the  differential  equations  (1)  and  boundary  conditions  (2-3).  To 
End  the  solution,  divide  the  time  interval  tf  —  to  into  M  subintervals  of  the  same  length 
h  =  (tf  —  to)/M.  Let  yk  =  y(to  +  kh )  be  the  value  of  the  solution  at  the  node  tk  =  t0  +  kh , 

k  —  0, . . . ,  M.  At  the  middle  point  between  tk- 1  and  tk,  denoted  by  tk- 1/2  =  tk  —  h/ 2,  the 

differential  equations  (1)  are  approximated  by  central  finite  difference: 

~  y^  1)  =  f  (tk~  1/2, Vk  2'  ^  (4) 

Or  equivalently,  yk  and  yk-\  are  constrained  by  the  equation 

Ekivki yk- 1)  =  yk  -  Vk-\  -  hf  (tk- 1/2, Vk  2k  i^)  =  M  =  !,  •  •  • , M  -  !•  (5) 

In  addition,  the  boundary  conditions  are  denoted  by 

E0(y0)  =  B0(y0 )  =  0 

em{vm )  =  B/{yf)  =  0 
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(6) 

(7) 


Treat  Y  =  (y$  y[  ■  ■  ■  Um)T  £  R2n(M+1)  as  the  unknowns.  The  same  number  of  equations 
are 

E(Y)  =  0  (8) 

where  E  =  ( E q  E\  .  .  .  E^)1'.  Now  the  problem  becomes  a  root-finding  problem  for  a  system 
of  2n(M  +  l)  nonlinear  algebraic  equations  (8).  ft  has  been  rigorously  established  that  under 
certain  conditions  on  smoothness  and  the  boundary  conditions,  the  following  holds  true20 

1.  Each  of  the  original  TPBVP  and  the  finite  difference  problem  has  a  unique  solution; 

2.  The  solution  of  the  above  finite  difference  problem  yk  is  a  second-order  approximation 
to  the  solution  of  the  TPBVP  y*(t)  at  each  tjtl  i.e., 

\\y*(tk)-yk\\  =  0(h2),  k  =  1,2, ...  ,M  (9) 

where  lim h,^0  0(h2) / h2  <  oo. 

For  ascent  guidance  applications,  since  the  time-to-go  tf  —  to  is  decreasing,  the  accuracy 
of  the  finite  difference  solution  will  be  higher  and  higher  as  h  becomes  smaller  even  for  a 
small  to  moderate  number  of  nodes. 

For  problems  with  free  final  time  tf ,  see  Section  5.1  for  a  treatment. 


5.2  Modified  Newton  Algorithm 


A  modified  Newton  method21  is  chosen  as  the  algorithm  for  solving  the  problem  (8)  for  its 
balance  between  complexity  and  effectiveness.  Starting  from  an  initial  guess  Y o  (which  may 
be  reliably  and  quickly  generated  from  vacuum  solution  using  the  algorithm  to  be  described 
in  Section  6),  the  search  direction  dj  in  the  j-th  iteration  is  determined  by  solving  the  linear 
algebraic  equations 

dj-^EiYj.t).  j  =  1,2,...  (10) 

Next,  determine  the  step  size  parameter  <jj  by  the  following  criterion 


dE{Yj- 1) 

BY 


co- 


max 

0<i 


(11) 


In  other  words,  starting  from  a3  =  1,  oq  is  halved  repeatedly  if  necessary  until  the  above 
condition  is  satisfied.  Then  the  update  is  given  by 


Yj  —  Y j_\  +  Cjdj,  0  <  Gj  <  1,  (12) 

This  choice  of  the  search  step  size  ensures  that  the  sequence  {||.E(V:;)||}  is  monotonically 
decreasing.  Convergence  is  achieved  when  \\E( Y 3 ) 1 1  is  no  greater  than  a  preselected  toler¬ 
ance.  The  possible  additional  function  evaluations  required  in  checking  the  above  step  size 
condition  pose  negligible  computational  burden  because  function  evaluations  are  not  expen¬ 
sive  in  this  setting.  The  result  on  the  other  hand  is  a  much  more  robust  algorithm,  especially 
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when  the  initial  guesses  are  not  close  to  the  final  solution.  The  step  size  selection  (11)  is  a 
deciding  factor  for  the  success  of  the  finite  difference  approach  in  solving  the  optimal  ascent 
problem. 

The  evaluation  of  the  Jacobian  dE/dY  can  certainly  be  done  analytically.  But  we 
believe  that  simple  numerical  finite  differencing  is  more  advantageous  in  this  case.  This  is 
because  again  unlike  in  the  cases  where  integrations  of  differential  equations  are  involved  for 
each  function  evaluation,  the  function  evaluations  here  are  purely  algebraic  and  fast.  Using 
analytical  Jacobian  offers  no  clear  computational  speed  benefits.  In  our  comparison  study 
we  have  seen  that  the  numerical  Jacobians  and  analytical  Jacobians  match  between  the  6th 
to  8th  digit.  With  the  scalings  described  in  this  report,  this  type  of  precision  appears  to 
be  more  than  adequate  for  convergence  to  occur.  On  the  other  hand,  analytical  Jacobian 
will  make  the  computer  code  significantly  more  complicated  because  second-order  partial 
derivatives  of  the  right  hand  sides  of  the  state  equations  are  needed.  Also,  when  some  of 
the  path  constraints  in  Eq.  (19)  become  active,  the  associated  Lagrange  multipliers  such 
as  the  one  in  Eq.  (39)  will  be  functions  of  state  and  costate,  adding  more  complexity  to 
the  analytical  Jacobian.  When  the  launch  vehicle  design  is  evolving  in  the  developmental 
stage,  or  if  the  same  ascent  guidance  code  is  desired  to  be  applied  to  different  launch  vehicle 
configurations,  fewer  labor-intensive  changes  to  the  code  would  be  required  with  numerical 
differencing.  One  exception  is  for  the  gradients  of  the  boundary  conditions  in  (7).  The 
gradients  for  these  conditions  are  evaluated  analytically  because  they  are  readily  available. 

At  first  glance,  solving  the  linear  system  (10)  may  seem  to  be  a  formidable  task  in  an 
on-board  environment  because  of  the  dimension  of  the  problem  (2n(M  + 1))  which  can  easily 
be  over  1000.  However,  a  closer  inspection  reveals  that  the  Jacobian  matrix  in  (10)  has  a 
special  sparse  pattern  clue  to  the  dependence  of  Ek  on  only  yk  and  yk_i  (cf.  Eq.  (5)). 
Therefore  a  very  efficient  algorithm,  both  in  speed  and  storage,  based  on  Gauss  eliminations 
and  sequential  back  substitutions  can  be  devised  to  solve  the  system  (10).  More  details  in 
this  regard  is  discussed  in  the  next  section. 


5.3  Gaussian  Elimination  and  Back  Substitution 

Expand  the  algebraic  system  (8)  in  first-order  Taylor  series  with  respect  to  small  changes 
AY” k  at  the  current  iterate  Y k.  At  an  interior  point,  k  —  1,  2, ...,  M,  this  gives 

Ek(Yk  +  ATfc,  Yk_ !  +  Ayw)  «  Ek(Yk ,  Yk_x)  +  — - —YYk_,  +  — ^A Yk  (13) 

OI  k- 1  O*  k 

For  a  solution  we  want  the  updated  value  E(Y  +  AY)  to  be  zero,  so  the  general  set  of 
equations  at  an  interior  point  can  be  written  as 

4^AYfc_!  +  ^A  Yk  =  - Ek(Yk ,  Yfc_0  (14) 

Oi  k- 1  C>¥  k 
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Similarly,  the  boundary  conditions  can  be  expanded  in  a  first-order  Taylor  series  for  incre¬ 
ments  that  improve  the  solution. 


d  E 
dY 

dEM 
dY  M 


-AY0 = - 

o 

■A  Ym  = 


E0{Y0) 

(15) 

Em(Y  m) 

(16) 

We  thus  have  in  Eqs  (14)-(16)  a  set  of  algebraic  linear  equations  to  be  solved  for  the 
corrections  A Y iterating  until  the  corrections  are  sufficiently  small.  The  equations  have  a 
special  structure,  because  the  Jacobians  involve  only  points  k  and  k  —  1.  Figure  5  illustrates 
a  typical  matrix  structure  of  the  complete  algebraic  equations  for  the  case  of  8  variables 
(two-dimensional  ascent  problem,  n  =  4)  and  4  mesh  points,  with  4  boundary  conditions 
imposed  at  each  of  the  endpoints.  In  the  figure,  “x”  represents  a  coefficient  of  the  FDEs, 
“v”  represents  a  component  of  the  unknown  solution  vector,  ,  “B”  is  a  component  of  the 
known  right-hand  side,  and  empty  spaces  represent  zeros. 
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Figure  5:  Matrix  structure  of  FDEs  with  boundary  conditions 


Because  of  the  special  “block  diagonal”  structure  of  the  matrix,  a  special  form  of  Gaus¬ 
sian  elimination,  which  can  minimize  storage  requirement  of  matrix  coefficients  by  packing 
the  elements  in  a  special  blocked  structure,  can  be  applied  to  solve  the  problem.22  General 
Gaussian  elimination  manipulates  the  algebraic  linear  equations  by  elementary  operations, 
such  as  dividing  rows  of  coefficients  by  a  common  factor  to  produce  unity  in  diagonal  ele¬ 
ments,  and  adding  appropriate  multiples  of  other  rows  to  produce  zeros  below  the  diagonal. 
In  this  special  Gaussian  elimination,  we  take  advantage  of  the  block  structure  by  performing 
a  bit  more  reduction  than  in  pure  Gaussian  elimination,  so  that  the  storage  of  coefficients 
is  minimized.  Only  a  small  subset  of  the  2 n(M  +  1)  x  2 n(M  +  1)  matrix  elements  needs  to 
be  stored  as  the  elimination  progresses.  Once  the  matrix  elements  reach  the  stage  in  Fig.  6, 
the  solution  can  follow  quickly  by  a  back  substitution  procedure. 
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Figure  6:  Target  structure  of  the  Gaussian  elimination 


The  entire  procedure  except  the  back  substitution  step,  operates  only  on  one  block  of 
the  matrix  at  a  time.  The  procedure  contains  4  types  of  operations: 

•  partial  reduction  to  zero  of  certain  elements  of  a  block  using  results  from  a  previous 
step; 

•  elimination  of  the  square  structure  of  the  remaining  block  elements  such  that  the  square 
section  contains  unity  along  the  diagonal,  and  zero  in  off-diagonal  elements; 

•  storage  of  the  remaining  nonzero  coefficients  for  use  in  later  steps; 

•  back  substitution. 


For  more  detail,  please  refer  to  Reference  22. 


Once  the  correction  AY  is  solved  from  the  algebraic  linear  equations,  we  update  the 
solution  by 

v  =  V,-.  +  °A  Yj  (17) 

Where,  the  subscript  “j”  denotes  the  j-th  iteration.  The  step  size  parameter  ay  is  determined 
by  the  following  criterion 


oy  =  max 

0<i 


1 

2* 


ET  (Y  7_i  + 


AY, 


)E(Y,  i  + 


AY, 


<  ET(YiAE(Yi_, 


(18) 
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Starting  from  =  1,  <jj  is  halved  repeatedly  if  necessary,  till  the  above  condition  is  sat¬ 
isfied.  This  choice  of  the  step  size  ensures  that  the  sequence  {||.E(1T;)||}  is  monotonically 
decreasing.  Convergence  is  achieved  when  ||.E(Tb)')||  is  no  greater  than  a  preselected  toler¬ 
ance.  This  step  size  selection  in  Eq.  (18)  is  a  deciding  factor  for  the  success  of  the  finite 
difference  approach  in  solving  the  optimal  ascent  problem.  The  possible  additional  function 
evaluations  required  in  checking  the  step  size  condition  pose  negligible  computational  burden 
since  function  evaluations  are  not  expensive  in  this  setting.  The  result  on  the  other  hand 
is  a  much  more  robust  algorithm,  especially  when  the  initial  guess  is  not  close  to  the  final 
solution. 


5.4  Continuation  on  Atmospheric  Density 

It  is  well  known  in  ascent  trajectory  optimization  that  the  strong  coupling  of  the  aerodynamic 
forces  with  the  orientation  of  the  body-axis  and  the  inequality  path  constraints  make  the 
convergence  of  any  algorithm  from  a  completely  “cold”  initial  start  difficult  to  achieve.  That 
is  the  reason  why  homotopic  filters,6  sometimes  up  to  4  levels  of  homotopy,7  are  used  to 
gradually  distort  the  solution  from  the  vacuum  solution  to  the  final  solution.  A  similar 
approach  is  taken  here,  except  that  the  finite  difference  approach  described  in  the  preceding 
section  appears  to  need  only  one  level  of  homotopy.  The  form  of  the  continuation  used  here 
is  simpler  and  is  only  on  the  atmospheric  density 

p  =  ep,  0  <  £  <  1  (19) 

For  each  fixed  £,  p  is  used  in  place  of  the  atmospheric  density  everywhere  p  appears  in  the 
state-,  costate-equations,  and  path  constraints.  Note  that  all  the  development  in  this  part 
of  the  report,  although  aimed  at  atmospheric  ascent,  applies  without  any  modifications  to 
vacuum  flight  when  p  =  0.  The  homotopic  parameter  £  is  initiated  at  0  for  the  vacuum 
solution,  and  gradually  increased  to  unity  for  full  atmospheric  solution.  The  converged 
solution  for  the  current  value  of  e  serves  as  the  starting  point  for  the  solution  with  the  next 
value  of  £  till  e  —  1.  We  have  been  using  an  increment  of  0.1  on  £,  which  seems  to  be 
adequate  in  our  cases.  In  on-board  applications,  no  continuation  is  usually  necessary  when 
the  solution  in  the  previous  guidance  cycle  is  used  as  the  starting  point  (that,  £  =  1). 

For  the  vacuum  solution,  in  many  cases  linear  interpolations  between  the  initial  values 
and  targeted  final  values  of  the  state  and  constant  guesses  for  pit)  suffice  to  be  an  initial 
guess  Y o  that  will  lead  to  convergence.  An  analytical  method  to  generate  optimal  vacuum 
ascent  trajectories  is  presented  in  Appendix  B.  This  method  is  based  on  a  succinct  summary 
by  Calise  et  a!6  of  a  number  of  elegant  results  on  optimal  vacuum  guidance  developed  in  the 
last  3  decades.  Some  modifications  have  been  made  in  the  version  in  Appendix  B.  With  this 
type  of  algorithm,  an  optimal  vacuum  solution  can  be  found  very  quickly  and  reliably.  This 
solution  can  in  turn  be  used  as  the  starting  solution  in  the  above  continuation  process. 
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Verification,  Validation  and  Testing 


A  number  of  test  cases  are  presented  in  this  section.  The  purpose  is  to  demonstrate  the 
working  of  the  endo-atmospheric  algorithm  developed  above.  More  extensive  testing  results 
are  contained  in  Part  IV  of  this  report.  All  the  cases  use  the  vehicle  parameters,  mass  prop¬ 
erty,  propulsion  system  Modeling,  and  aerodynamic  Modeling  of  the  X-33  vehicle.  The  X-33 
is  a  single-stage  suborbital  vehicle  with  a  lifting-body  configuration  (see  Fig.  3).  To  make 
the  vehicle  orbit-capable,  the  specific  impulse  is  doubled  in  the  simulations.  All  the  missions 
are  launched  from  the  Kennedy  Space  Center  (KSC).  The  open-loop  solutions  obtained  us¬ 
ing  the  method  described  in  this  report  are  compared  with  those  obtained  by  using  other 
methods.  Closed-loop  simulations  are  performed  to  assess  the  feasibility  of  on-board  ascent 
guidance  with  the  approach  proposed  in  this  report.  Unless  stated  otherwise,  all  the  vectors 
are  expressed  in  the  inertial  launch  plumbline  coordinate  system  introduced  in  Section  2.2. 


6.1  Terminal  Conditions  and  Final  Time  Adjustment 


The  finite  difference  approach  described  in  Section  4.1  is  most  convenient  for  fixed  final-time 
problems.  Similar  to  what  is  done  in  Ref.  [7],  the  optimal  ascent  problem  is  solved  as  a 
series  of  fixed  final-time  problems  to  maximize  the  orbital  energy,  i.e.,  to  minimize 


J  = 


1 

rf 


(1) 


The  final  time  is  adjusted  sequentially  until  the  optimal  value  of  J  is  equal  to  the  specified 
(negative)  orbital  energy. 

Four-Constraint  Problem 

The  orbital  insertion  conditions  are  given  by  the  radius  r  j,  velocity  VJ ,  orbital  inclination 
i*,  and  flight  path  angle  7 j.  Note  that  7^  need  not  necessarily  be  zero  in  this  case.  These 
conditions  are  equivalent  to  specified  semi-major  axis,  eccentricity,  orbital  inclination,  and 
true  anomaly  at  insertion  point.  Let  ljv  be  a  unit  vector  parallel  to  the  polar  axis  of  the 
Earth  and  pointing  to  the  North.  For  each  fixed  tf  the  3  orbital  insertion  conditions  (15)  in 
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this  case  can  be  expressed  as 

\rT,r  f~\r?  =  0  (2) 

1w'(r/  x  V f)  —  \\rf  x  Vf\\  cosi*  =  0  (3) 

rTfVf  ~rfvfs[nlf  =  0  (4) 

For  simplicity,  we  replace  ry  in  (1)  with  r*j.  Taking  dot  products  of  the  transversality 
condition  Eq.  (25)  with  Vf,  Eq.  (26)  with  r/,  Eq.  (25)  with  V /,  and  both  Eqs.  (25)  and 
(26)  with  hf  =  r f  x  V f,  we  can  eliminate  the  multiplier  vector  v  and  obtain 

+  =  0  (5) 

VTSPv,~Vf  =  0  (6) 

{hTfPrf)[hTf(rf  x  ljv)]  +  (hTfpVf)[hTf(V f  x  1N)]  =  0  (7) 

The  above  6  equations  (2-7)  constitute  the  terminal  boundary  conditions  (3)  for  this  case. 

To  adjust  the  final  time,  we  have  found  that  the  secant  method  works  very  well  for 
this  purpose.  Let  ^  and  t^"1  be  two  consecutive  estimates  of  tf  used  to  solve  the  above 

problem,  and  and  be  the  values  of  J  when  and  are  used  to  solve  the 

optimal  ascent  problem.  Then  the  next  choice  of  tf  is  given  by 

d+I1  =  C  -  (jtlll-l))  (jm  -  •n,  k  =  1,2, ...  (8) 


The  correct  tf  is  found  when  —  J*  |  is  within  a  preset  tolerance,  and  the  final  velocity 
Vf  will  be  within  a  small  neighborhood  of  Vf.  To  use  the  secant  search  (8),  two  starting 

values  of  tj1  and  tj  }  are  needed.  For  on-board  guidance,  the  converged  value  of  tf  in  the 

previous  guidance  cycle  is  the  logical  choice  of  t ^ .  Since  the  guidance  solution  in  the  current 
cycle  will  be  very  close  to  the  previous  one,  the  variational  relationship  of  5J  =  H(tf)8tf  in 
optimal  control  problems  can  be  used  to  generate  tj' 


Ml  =  JO)  J(01  -  £ 

’  ’  H(tf) 


where  H(t is  the  final  value  of  the  Hamiltonian.  It  is  possible  to  simply  use  (9)  to  search 

for  all  tj  ^ ,  k  >  1,  as  in  Refs.  [7-8].  But  we  have  found  the  secant  method  (8)  to  be  much 
more  robust  and  efficient  than  the  exclusive  use  of  Eq.  (9). 


Five- Constraint  Problem 

The  orbital  insertion  conditions  are  given  by  the  radius  r*f,  velocity  Vf,  orbital  inclination 
i* ,  flight  path  angle  7])  =  0  (insertion  at  perigee/ apogee,  or  into  a  circular  orbit),  and 
longitude  of  ascending  node  fT.  A  mission  to  rendezvous  with  an  orbiting  spacecraft  would 
call  for  the  5th  constraint  on  f2*.  With  i*  and  specified,  the  direction  of  the  angular 
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momentum  vector  of  the  orbit  is  fixed,  given  in  an  Earth  centered  inertial  (ECI)  system  by 
the  unit  vector  lfcy  =  (sinf2*  sin i*,  —  cosf2*  sin i*,  cosi*)T.  Let  1  ^  be  the  representation  of 
if67  in  the  inertial  launch  plumbline  system.  The  4  orbital  insertion  conditions  for  a  fixed 
final  time  are 


1  *2 
~2rf 

=  0 

(10) 

rJVf 

=  0 

(11) 

rf ^ 

=  0 

(12) 

VTflh 

=  0 

(13) 

By  taking  dot  products  of  the  transversality  condition  Eq.  (25)  with  V f,  Eq.  (26)  with 
V  f,  then  with  rf,  and  making  use  of  the  above  conditions,  we  again  eliminate  the  multiplier 
vector  v  and  obtain 


(yTfprf)r)  -  {rTfPvf)Vj  =  0  (14) 

VTfPVf-Vf  =  0  (15) 

The  final  time  tf  is  adjusted  the  same  way  by  Eq.  (8)  so  that  Vf  =  VJ . 

For  staged  launch  vehicles  where  all  the  stages  except  for  the  last  one  burn  to  the  de¬ 
pletion  of  the  fuel,  the  burn  time  of  each  stage  is  known  for  given  throttle.  The  approach 
in  this  report  is  still  applicable  with  minor  modifications.  Note  that  the  break  points  in  the 
finite  difference  discretization  do  not  have  to  be  equally  spaced.  In  the  case  of  a  IV-stage 
vehicle,  N  different  values  for  the  step  parameter  h  may  be  selected,  each  for  a  stage,  such 
that  the  staging  times  are  located  exactly  at  the  last  nodes  of  the  stages.  The  changes  of 
mass  and  thrust  due  to  staging  can  then  be  accommodated  easily.  The  burn  time  of  the  last 
stage  is  adjusted  to  meet  the  final  velocity  condition. 

A  number  of  other  orbital  insertion  conditions  were  also  examined  and  implemented.  But 
the  principle  is  the  same,  thus  we  will  not  repeat  them  in  this  report. 


6.2  Open-loop  Solutions 

The  open-loop  solutions  are  generated  using  the  finite-difference  (FD)  method.  For  verifica¬ 
tion  and  validation  of  the  approach,  the  trajectories  are  compared  to  the  ones  obtained  with 
other  established  trajectory  optimization  tools.  The  wind  velocity  Vw  is  set  to  be  zero  in 
the  solutions.  The  target  orbit  is  a  circular  orbit  at  the  altitude  of  185.2  km  (100  nm).  Two 
orbital  inclinations  are  used:  i*  =  51.6  deg  (the  orbital  inclination  of  the  International  Space 
Station),  and  i*  =  28.5  deg  (the  minimum  orbital  inclination  from  KSC).  A  third  case  is  the 
same  circular  orbit  with  i*  =  51.6  deg  and  Q*  =  —104  deg  (a  5-constraint  problem).  The 
initial  conditions  correspond  to  those  after  5  seconds  of  vertical  ascent  (to  clear  the  tower) 
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with  a  given  launch  mass  of  the  X-33.  The  following  path  constraints  are  imposed 


(16) 

(17) 


\qa\  <  626.74  N-rad/m2  (750  psf-deg) 

T  <  4.0  (g) 

The  engine  throttle  is  set  at  unity  before  the  thrust  acceleration  constraint  is  active.  The 
constraint  on  dynamic  pressure  will  be  added  later  in  closed-loop  simulations  using  the  tech¬ 
nique  presented  in  Section  3.7.  In  the  following  results,  100  nodes  were  used  in  FD  solutions. 
For  comparison,  the  solutions  for  the  same  missions  were  also  obtained  by  a  collocation  soft¬ 
ware,  Sparse  Optimal  Control  Software  (SOCS),1'  and  by  a  trajectory  optimization  software 
based  on  a  pseudo-spectral  method,  Direct  and  Indirect  Dynamic  Optimization  (DIDO).18 
The  SOCS  solutions  had  the  same  100  nodes,  and  the  DIDO  solutions  had  20  nodes.  The 
performance  index  for  both  SOCS  and  DIDO  was  the  final  time. 

Table  1  summarizes  the  flight  times  tf  and  final  mass  rrif  of  the  solutions.  The  results 
for  all  the  missions  under  different  methods  are  very  close.  The  slight  differences  among  the 
solutions  were  mostly  due  to  the  discretization  errors  of  the  different  schemes  used  in  the  3 
methods.  This  comparison  clearly  supports  the  validity  of  the  FD  approach.  The  difference, 
however,  is  in  computation  time.  Depending  on  initial  guesses,  SOCS  and  DIDO  could  use 
tens  of  minutes  to  find  the  solution,  whereas  the  FD  method  takes  only  a  small  fraction  of 
one  second  to  find  the  solution. 

Figure  7  shows  the  3-D  ascent  trajectories  and  ground  tracks  for  the  cases  of  i*  =  51.6 
deg  and  0  =  —104  deg  in  the  inertial  launch  plumbline  coordinate  frame.  As  expected,  the 
trajectory  with  ascending  node  constraint  had  larger  out-of-plane  motion.  Figure  8  depicts 
the  angle  of  attack,  pitch  angle  and  yaw  angle  along  all  the  3  trajectories  by  the  FD  ap¬ 
proach.  The  constraint  on  the  ascending  node  usually  is  one  that  requires  larger  yaw  and 
roll  maneuvers  and  thus  is  a  demanding  constraint.  Indeed  larger  yaw  angle  is  observed  in 
Fig.  8  for  that  case.  The  variations  of  qa  and  axial  thrust  acceleration  are  plotted  in  Fig.  9 
Both  path  constraints  were  accurately  enforced. 


6.3  Closed-loop  Simulations 

In  closed-loop  simulations,  the  point-mass  vehicle  dynamics,  atmospheric  Modeling,  propul¬ 
sive  and  aerodynamic  forces,  and  winds  were  simulated  by  a  FORTRAN  program.  The  FD 
algorithm,  which  is  also  implemented  in  FORTRAN,  was  called  once  every  second  to  recal¬ 
culate  the  optimal  solution  based  on  the  current  condition.  The  trajectory  simulations  used 
the  first  data  in  the  optimal  body-axes  attitude  and  engine  throttle  solutions  (corresponding 
to  current  time)  as  the  guidance  commands.  No  delays  or  actuator  dynamics  were  simulated. 
The  missions  and  initial  conditions  were  the  same  as  in  the  open-loop  solutions.  In  addition 
to  path  constraints  (16)  and  (17),  a  dynamic  pressure  constraint  is  also  imposed 

q  <  18194.4  N/m2  (380  psf)  (18) 

Closed-loop  Simulations  without  Winds 
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The  previous  3  missions  were  repeated  in  closed-loop  simulations  with  the  wind  velocity 
Vw  set  to  be  zero.  The  dynamic  pressure  constraint  was  handled  using  the  technique  de¬ 
scribed  in  this  report.  Thirty  nodes  were  used  in  the  FD  guidance  solution.  For  on-board 
guidance,  fewer  nodes  can  be  used  than  in  an  open-loop  solution  because  the  time-to-go  is 
decreasing,  thus  the  accuracy  of  the  solution  is  increasing  while  the  computation  demands 
would  not  be  unnecessarily  high.  A  converged  solution  was  loaded  before  “launch”  as  the 
starting  point.  On  a  desktop  computer  with  a  1  GHz  processor,  the  CPU  time  that  each 
guidance  solution  call  took  ranged  from  0.07  to  0.25  second,  without  any  optimization  or 
streamlining  of  the  code.  All  the  orbital  insertion  conditions  were  met  accurately  (See  Table 
3  in  the  following).  Table  2  lists  the  flight  times  and  final  masses  for  all  the  3  missions.  Also 
in  the  table  are  open-loop  solutions  by  SOCS  with  the  addition  of  the  dynamic  pressure 
constraint  (18).  In  all  3  cases  the  closed- loop  guided  trajectories  and  the  SOCS  open-loop 
solutions  have  differences  in  the  final  mass  of  about  120  kg  or  less.  These  small  discrepan¬ 
cies  are  largely  attributed  to  the  differences  between  closed-loop  simulations  and  open-loop 
trajectories. 

Figure  10  shows  the  altitude  versus  inertial  velocity  along  the  three  trajectories.  The 
variations  of  dynamic  pressure,  throttle,  and  axial  thrust  acceleration  are  plotted  in  Fig.  11. 
The  effectiveness  of  the  feedback  approach  described  in  Section  3.7  for  enforcing  the  dynamic 
pressure  constraint  is  clearly  seen:  the  dynamic  pressure  constraint  (18)  is  accurately  met 
in  all  cases  and  the  feedback  throttle  modulations  are  smooth.  This  is  much  similar  to  the 
” throttle  bucket”  the  Shuttle  employs  for  the  same  purpose,  except  that  the  Shuttle  throttle 
bucket  is  pre-programmed  (open  loop)  prior  to  launch.  In  comparison,  the  SOCS  solutions 
all  have  considerable  zigzags  in  the  throttle  when  the  constraint  (18)  is  active  (not  shown 
in  the  figure).  From  Table  2,  the  performance  (in  terms  of  final  mass)  of  the  closed-loop 
trajectory  with  the  feedback  adjustment  of  the  throttle  is  apparently  about  the  same  as  that 
of  the  open-loop  optimal  solution.  The  variations  of  qa  and  a  are  in  Fig.  12.  Figure  13 
contains  the  histories  of  pitch,  yaw  and  roll  angles  along  the  3  trajectories.  As  expected, 
the  mission  with  ascending  node  constraint  called  for  more  yaw  and  roll  maneuvers  than  the 
other  two  missions. 

Closed-loop  Simulations  with  Winds 

So  far  in  the  simulations  the  wind  velocity  has  been  assumed  to  be  zero.  To  test  the 
ascent  guidance  algorithm  in  a  more  realistic  setting,  winds  were  added  in  the  simulations. 
The  wind  profiles  were  based  on  the  measured  wind  velocities  at  different  altitudes  at  KSC, 
and  then  smoothed  for  guidance  purposes.  In  each  simulated  trajectory,  the  same  smoothed 
wind  profile  was  used  in  both  the  guidance  solution  and  simulation  of  the  launch  vehicle 
dynamics.  We  realize  that  the  actual  wind  the  launch  vehicle  experiences  will  likely  be 
different  from  the  measured  wind  because  of  the  time  delay  between  the  launch  and  the  time 
when  the  measurement  was  taken.  But  the  purpose  here  is  to  demonstrate  how  the  ascent 
guidance  algorithm  would  perform  in  the  presence  of  winds  should  perfect  wind  information 
be  available.  We  have  also  tested  the  cases  when  the  wind  profiles  used  in  guidance  algorithm 
were  correlated  but  not  exactly  the  same  as  the  wind  profiles  in  simulations.  In  those  cases 
the  enforcement  of  path  constraints  is  affected  by  how  well  the  measured  wind  profiles  match 
the  actual  wind  profiles.  With  the  volatile  winds,  some  of  which  are  quite  strong  (up  to  75 
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m/s),  the  limit  of  qa  needs  to  be  higher  in  order  to  fly  through  the  winds.  Thus  the  constraint 
(16)  is  changed  to 

\qa\  <  1671.32  N-rad/m2  (2000  psf-deg)  (19) 

Ten  trajectories  were  simulated  for  each  of  the  above  three  missions,  with  a  different  wind 
profile  applied  in  each  trajectory.  Unlike  open-loop  ascent  guidance,  no  “pre-launch”  ad¬ 
justments  or  updates  were  required  for  the  guidance  algorithm  for  each  simulation.  The 
only  update  for  each  flight  was  the  wind  profile.  The  closed-loop  ascent  guidance  algorithm 
automatically  incorporates  the  wind  data  in  the  guidance  solution.  The  starting  point  was 
still  the  same  zero-wind  solution  used  before.  The  CPU  time  for  each  guidance  cycle  did  not 
differ  from  the  previously  recorded  value  by  any  noticeable  margin. 

The  simulation  results  are  summarized  in  Table  3.  The  first  6  quantities  with  A  are  the 
orbital  insertion  condition  errors.  The  minimum,  average,  and  maximum  values  of  these 
errors  among  the  10  trajectories  for  each  mission  are  listed.  The  altitude  errors  were  all 
within  0.2  meter  and  velocity  errors  within  0.4  m/s.  The  quantity  A  a  is  the  error  in  semi¬ 
major  axis  of  the  final  orbit,  which  is  a  parameter  combining  the  effects  of  A r/  and  A Vf. 
The  maximum  of  A  a  was  no  greater  than  0.6  km.  All  other  errors  were  very  small  as  well. 
The  peak  values  of  \qa\  were  essentially  all  within  the  specified  range.  The  product  of  q/3 , 
where  f3  is  the  sideslip  angle,  was  not  listed  for  it  remained  practically  zero  because  of  the  roll 
maneuvers  by  the  launch  vehicle  to  “fly  into  the  wind” .  Even  when  the  wind  profiles  used 
in  simulations  were  different  from  but  correlated  to  the  wind  profiles  used  in  the  guidance 
solution,  such  as  in  actual  launch,  the  peak  values  of  \q/3\  in  general  were  still  significantly 
smaller  than  those  of  \qa\  (not  shown  here).  The  average  final  masses  for  all  the  3  missions 
were  quite  close  to  those  in  Table  2  for  the  trajectories  without  winds.  It  should  be  mentioned 
that  whether  or  not  the  wind  data  is  included  in  the  optimal  guidance  solutions  can  make 
a  sizable  difference  in  the  performance.  If  the  wind  profiles  are  merely  used  outside  the 
guidance  solution  to  limit  the  attitude  guidance  commands  for  observing  the  ga-constraint, 
the  results  using  the  X-33  vehicle  Model  can  mean  up  to  500  kg  less  payload  delivered  into 
the  same  orbit. 


Table  1:  Open- loop  performance  comparison 


Mission 

FD 

SOCS 

DIDO 

tf  (sec) 

mf  (kg) 

tf  (sec) 

mf  (kg) 

tf  (sec) 

771/  (kg) 

i  =  51.6° 

320.48 

38627.5 

321.00 

38476.2 

321.12 

38471.8 

i  =  28.5° 

317.73 

39135.7 

318.05 

39011.0 

318.18 

39007.2 

O  =  -104° 

322.71 

38220.1 

322.33 

38236.5 

322.81 

38166.9 
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Table  2:  Closed-loop  simulations  without  winds 


Mission 

FD 

SOCS  (open-loop) 

tf  (sec) 

mf  (kg) 

tf  (sec) 

mf  (kg) 

i  =  51.6° 

326.57 

38182.4 

323.23 

38293.6 

i  =  28.5° 

323.59 

38706.7 

320.24 

38828.6 

0  =  -104° 

327.33 

38023.9 

324.5 

38054.6 

Table  3:  Closed-loop  simulations  with  winds  (10  wind  profiles  for  each  mission) 


i  =  51.6° 

i  =  28.5° 

n  =  -104° 

Min 

Mean 

Max 

Min 

Mean 

Max 

Min 

Mean 

Max 

A rf  (m) 

-0.0165 

0.0128 

o.ior 

-0.0115 

0.0152 

0.112 

-0.0145 

0.0484 

0.117 

A Vf  (m/s) 

0.1531 

0.2290 

0.3007 

0.1625 

0.2352 

0.3327 

0.1583 

0.1909 

0.25388 

A  a  (km) 

0.2581 

0.3856 

0.5063 

0.2735 

0.3961 

0.5602 

0.2665 

0.3215 

0.4275 

A i  (deg) 

2.4E-5 

4.6E-5 

5.9E-5 

-1.9E-06 

-4.2E-07 

5.7E-07 

-5.8E-4 

-2.4E-4 

7.9E-05 

AS2  (deg) 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

-7.2E-4 

-2.9E-4 

9.8E-05 

Ay  (deg) 

-3.3E-4 

1.2E-4 

1.7E-3 

-2.6E-4 

1.3E-4 

1.5E-3 

-3.0E-4 

6.4E-4 

1.7E-3 

max  \qa\* 

919.4 

1268.3 

1608.3 

694.2 

1157.8 

1641.8 

1085.3 

1255.8 

1690.3 

mf  (kg) 

38010.4 

38131.1 

38238.5 

38507.2 

38626.3 

38785.4 

37789.4 

37948.9 

38065.3 

tf  (sec) 

326.51 

328.66 

331.21 

322.68 

326.77 

329.51 

327.29 

329.66 

333.19 

*  unit  in  N-rad/m2 


inclination  =  516.deg 


Figure  7:  Three-dimensional  ascent  trajectories  in  launch  plumbline  system 
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Figure  8:  Angle  of  attack,  pitch  and  yaw  angles  in  open-loop  solutions 
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Figure  9:  Variations  of  qa  and  axial  thrust  acceleration  in  open-loop  solutions 


Figure  10:  Altitude  versus  inertial  velocity  along  closed-loop  trajectories 
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Figure  11:  Dynamic  pressure,  throttle  and  axial  acceleration  along  closed-loop  trajectories 
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Figure  12:  Variations  of  qa  and  a  along  closed-loop  trajectories  (no  winds) 
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Figure  13:  Variations  of  Euler  angles  along  closed-loop  trajectories  (no  winds) 
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7 


Exo- Atmospheric  Closed-Loop  Ascent 
Guidance 

7.1  Introduction 

Launch  ascent  mission  planning  and  guidance  is  an  engineering  area  where  routine  appli¬ 
cations  of  optimization  tools  and  optimal  control  theory  are  a  necessity  rather  than  nicety. 
Indeed,  on-board  algorithms  for  solving  the  optimal  ascent  problem  are  the  foundation  of 
closed-loop  ascent  guidance  for  upper  stages  of  launch  vehicles  since  the  1960s.  The  exam¬ 
ples  of  classical  exo-atmospheric  optimal  ascent  guidance  algorithms  include  the  Iterative 
Guidance  Mode  (IGM)  guidance  employed  for  the  Saturn  V  rockets1  and  the  Powered  Ex¬ 
plicit  Guidance  (PEG)  for  the  Space  Shuttle.2  When  the  optimal  ascent  problem  is  solved 
on-line  repeatedly  with  the  current  condition  as  the  initial  condition,  the  guidance  solution 
is  in  effect  closed  loop. 

For  multiple-stage  launch  vehicles,  it  is  well  known  that  allowing  a  coast  arc  of  optimal 
duration  between  two  upper  stages  can  reduce  propellant  consumption  (or  increase  orbital 
insertion  mass),  in  some  cases  quite  appreciably.  In  other  scenarios  a  long  coast  is  necessary 
between  two  burns,  such  as  insertion  into  high-altitude  orbits,  or  orbital  transfers.  However, 
there  are  additional  challenges  resulting  from  the  presence  of  coast  arcs  and  multiple  burns. 
Multiple  burns  and  coast  arcs  typically  render  the  optimal  control  problem  much  more  sen¬ 
sitive  and  increase  the  difficulty  in  achieving  convergence  of  the  algorithm.  This  heightened 
sensitivity  is  largely  clue  to  two  reasons:  (a)  long  coast  arcs  tend  to  amplify  any  variations 
in  the  condition  at  the  beginning  of  the  coast;  (b)  the  final  stage/burn  usually  has  a  much 
smaller  thrust  acceleration,  and  it  would  be  incapable  of  making  large  flight  path  corrections 
(indeed,  the  last  burn  in  many  orbital  insertion  cases  is  almost  tangential,  mainly  responsi¬ 
ble  for  pushing  the  velocity  to  the  required  value).  The  fact  that  even  today,  missions  with 
multiple-burns  and  coast  arcs  are  planned  and  executed  by  the  ground  control  is  probably 
attributable  in  no  small  part  to  these  challenges.  The  recent  on-going  efforts  in  striving 
to  achieve  the  capabilities  of  responsive  space  launch  mission  planning  and  operations  and 
autonomous  space  operations  have  brought  renewed  strong  interest  in  this  subject. 

To  overcome  the  aforementioned  difficulties,  researchers  have  continued  to  seek  more 
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robust  and  reliable  methods  to  handle  the  optimal  ascent  problem  with  multiple  burns  and 
coast  arcs.  References6, 10,11,23-26  represent  only  a  small  sampling  of  the  literature.  With  the 
exception  of  the  work  by  Dukeman, 10,11  most  existing  work  treats  the  problem  essentially 
in  a  single-shoot  formulation,  i.e.,  the  missing  initial  costate  is  iterated  on  to  satisfy  the 
terminal  constraints  and  transversality  conditions.  References10  and11  employ  a  numerical 
multiple-shooting  approach  to  the  problem  including  the  endo-atmospheric  ascent  portion. 

In  this  report  we  present  an  analytical  multiple-shooting  method  for  exo-atmospheric 
multi-burn  optimal  ascent.  The  overarching  theme  in  this  effort  is  to  strive  for  enhanced  ro¬ 
bustness  and  reliability  of  the  algorithm.  In  this  approach  the  optimal  trajectory  is  treated 
in  segments  of  burn  and  coast  arcs.  The  costate  in  each  segment  is  expressed  in  closed- 
form  solution  and  the  state  in  analytical  solution  involving  thrust  quadratures.  Such  a 
multiple-shooting  formulation  reduces  the  sensitivity  of  the  problem  with  respect  to  the  pa¬ 
rameters  to  be  found,  especially  beneficial  to  the  cases  with  long  coast  arcs.  The  orbital 
insertion  conditions,  transversality  conditions,  and  matching  continuity  conditions  at  the 
break  points  constitute  a  system  of  nonlinear  algebraic  equations  with  an  analytical  Jaco¬ 
bian.  For  numerical  solutions,  the  existing  work  almost  invariably  uses  the  classical  or  some 
modified  version  of  the  Newton  method.  Here  we  adopt  the  highly-regarded  Powell’s  dog-leg 
method28  for  improved  robustness  and  reliability  of  the  overall  algorithm  in  more  difficult 
and  even  possibly  singular  cases  (where  the  optimal  coast  time  reduces  to  zero).  Another 
major  component  of  this  report  is  a  detailed  analysis  of  the  class  of  optimal  ascent  problems 
under  discussion.  This  analysis  reveals  several  properties  previously  unknown  to  the  commu¬ 
nity,  to  the  best  of  our  knowledge.  The  results  allow  us  to  resolve  satisfactorily  a  numerical 
scaling  mismatch  issue  in  the  optimal  ascent  problem  that  could  otherwise  severely  hinder 
convergence  reliability  of  the  algorithm.  All  these  measures  contribute  to  a  robust  and  fast 
algorithm  for  optimal  multi-burn  ascent  planning  and  closed-loop  guidance.  The  algorithm 
is  verified  with  an  industry  standard  software,  Optimal  Trajectories  by  Implicit  Simulation 
(OTIS).29 

7.2  Multi-Stage  Optimal  Ascent  Problem  with  Coast 

Instead  of  seeking  maximum  generality,  we  choose  tacitly  to  present  the  development  in 
a  burn-coast-burn  pattern  for  readability  of  the  presentation.  There  is  no  methodological 
difficulty  to  extend  the  development  here  to  any  combinations  of  other  numbers  of  burn  and 
coast  arcs.  Borrowing  the  Shuttle  terminology,  we  will  refer  the  burnout  of  the  first  vacuum 
burn/stage  as  Main  Engine  Cut-Off  (MECO),  and  the  second  burn/stage  after  a  coast  as 
the  Orbital  Maneuver  System  (OMS)  burn.  For  on-orbit  maneuvers  it  is  understood  that 
both  burns  can  be  made  by  the  same  engine,  and  there  is  no  staging  at  the  end  of  the  first 
burn. 
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7.3  Problem  Formulation 


After  the  launch  vehicle  clears  the  dense  atmosphere,  the  three-dimensional  point-mass  equa¬ 
tions  of  motion  in  vacuum  are 


r  =  V 

(i) 

V  =  §0)  + 

m(t) 

(2) 

T 

rh  = - — 

9o^sp 

(3) 

where  r  G  R3  and  V  G  R3  are  the  position  and  inertial  velocity  vector  in  an  inertial  frame. 
The  gravitational  acceleration  g  is  a  function  of  r,  and  g0  is  the  magnitude  of  g  at  a  reference 
radius  Rq.  The  engine  thrust  magnitude  is  T,  and  the  unit  vector  1  t  defines  the  direction  of 
the  thrust  vector.  The  vehicle  mass  rate  m  is  determined  by  the  last  equation  above  where 
Isp  is  the  specific  impulse  of  the  engine.  An  ingenious  approximation  to  the  gravitational 
acceleration  is  the  so-called  linear  gravity  approximation30 
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'  -  2 

—  =  —  to  r 

r 


(4) 


where  g  is  the  gravitational  parameter  of  the  Earth,  f  is  another  reference  radius  (e.g.,  an 
average  value  of  r  along  the  ascent  trajectory),  and  Co  =  \J g/f3  is  the  Schuler  frequency  at 
f.  This  approximation  preserves  the  more  important  characteristic  of  the  gravity  in  ascent 
flight,  the  direction,  and  enables  analytical  solution  to  the  co-state  equation  in  optimal  exo- 
atmospheric  flight  as  shall  be  seen  later.  The  minor  gravity  magnitude  difference  caused  by 
this  approximation  in  ascent  typically  has  little  influence  on  the  validity  of  the  solution.  For 
multiple-burn  trajectory,  different  value  of  the  Co  may  be  used  for  each  burn  arc  so  that  any 
approximation  effect  on  gravity  magnitude  is  minimized. 

An  assumption  on  the  class  of  problems  treated  is  that  the  initial  condition  (r0,  Vo) 
is  known.  Another  practical  assumption  on  the  methodology  developed  in  this  report  is 
that  the  burn  times  of  all  the  powered  stages  except  for  the  last  are  all  determined  by  their 
propellant  loading  and  mass  flow  rate.  The  burn  time  of  the  last  powered  stage  and  the 
coast  times  are  optimized.  Therefore  the  final  time  for  the  complete  trajectory  is  free.  For 
two-burn  orbital  transfers,  the  burn  time  of  the  first  burn  may  also  need  be  optimized.  In 
such  a  case,  all  the  developments  and  equations  in  the  rest  of  this  report  remain  valid,  and 
one  additional  interior  switching  condition  will  be  added  (cf.  Ref.10). 

For  better  numerical  conditioning,  the  distances  are  normalized  by  i?0 >  the  velocities  by 
RRogoi  and  the  time  by  y/R0/g0.  With  some  abuse  of  notation,  we  will  still  use  r  and  V  to 
denote  the  dimensionless  position  and  velocity  vector,  respectively.  With  the  above  linear 
gravity  approximation,  the  dimensionless  equations  of  motion  are 


r'  =  V  (5) 

V1  =  —to2r  +  At1t  (6) 
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where  in  above  equations  the  differentiation  is  with  respect  to  the  nondimensional  time 
T  =  t/y/Ro/go,  LU  =  vVW)  3  is  the  nondimensional  Schuler  frequency,  and  At  =  T / mg o 
is  the  instantaneous  thrust  acceleration  in  g.  The  mass  rate  equation  becomes 

,  T 

rn  — - 

c 

where  c  =  Isp/ y/Ro/g f  is  a  constant  for  each  powered  stage.  The  orbital  insertion  conditions 
are  generally  in  terms  of  k  equality  constraints  (k  <  6)  on  the  final  state 


Hrf,Vf)  =  0 


(8) 


The  optimal  thrust  direction  vector  1  t  is  determined  by  the  solution  of  the  optimal 
control  problem  that  minimizes  the  performance  index 


J  = 


rTf 


mdiT  = 


rTf  rp 


dr 


(9) 


'  TO 


'  TO 


It  is  clear  that  this  performance  index  amounts  to  minimization  of  propellant-consumption 
for  a  given  initial  mass.  The  standard  optimal  control  theory14  calls  for  the  use  of  the 
Hamiltonian 


H  =  pJV  -  u2pyr  +  PylTAT  -  pm—  -  — 

c  c 

=  pjY-u2plr  +  T^P^lT  P 


mg  o 


—  -  -  :=  H0  +  TS 

c  c 


(10) 


In  above  Hamiltonian  pr  and  pv  constitute  the  costate  vector,  satisfying  the  differential 
equation 


pr 

Pv 


(  dHf  dr  \ 
\dH/dV  ) 


_  (u  Pv 

~Pr 


(11) 

KdH/ov ;  '  "  ' 

In  particular,  pv  is  called  the  primer  vector  because  the  optimal  thrust  direction  1  t  = 
Pv/WPvW-31  The  switching  function  S,  dehned  as 


5  = 


P^At  Pr, 


mg  o  c  c 

determines  when  to  use  full  thrust  and  when  to  coast  (T  =  0).  Specifically, 

T 

J.  yy 


(12) 


T  = 


max  if  ‘S'  >  0, 

0  if  S  <  0. 


(13) 


The  case  of  singular  thrust  arcs  (when  S'  =  0  in  a  finite  period  of  time  and  the  thrust  takes 
some  intermediate  values  0  <  T  <  Tmax )  is  not  considered  because  it  is  well  known  that, 
except  for  a  few  pathological  cases,  the  singular  solutions  are  not  optimal  in  exo-atmospheric 
flight.32  Suppose  in  our  discussion  that  the  coast  between  the  first  and  second  stage  ends 
at  a  to-be-determined  time  tqms ■  It  is  then  necessary  that  S(tqms )  =  0  (and  S'(r)  >  0 
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for  t  G  (■ tomSit, /))■  Since  the  total  final  time  T/  is  free  for  this  optimal  ascent  problem, 
the  Hamiltonian  H  defined  in  Eq.  (10)  is  equal  to  zero  along  the  optimal  trajectory.14  As  a 
result,  the  condition  S(toms)  =  0  is  equivalent  to 

Hq{toms )  =  0  (14) 

This  equation  avoids  the  computation  of  the  mass  costate  pm  thus  is  preferred  to  the  condition 
S(toms )  =  0.  For  multiple  coast  arcs,  refer  to  Ref.10  for  additional  interior  point  conditions. 
In  light  of  this  observation  and  other  necessary  conditions  discussed  above,  we  conclude  that 
the  mass  costate  pm  need  not  be  explicitly  computed. 


7.4  Analytical  Solution  for  Burn  Arcs 


The  analytical  solutions  to  the  costate  equation  (11)  and  analytical  solution  to  state  equa¬ 
tions  (5-6),  first  presented  in  Ref.30  and  later  also  applied  in  Refs.5  and,12  are  given  below 
for  completeness.  Suppose  that  the  starting  time  for  the  ascent  trajectory  is  r0.  Let  pVo  and 
prQ  are  the  (to-be-determined)  initial  conditions  for  the  costate  at  tq.  Define 


A(r) 


Pv(r)  \ 

-pr{r)/uj  )  ’ 


A0  — 


PV0  \ 

~PrJu  J 


For  t  >  To  the  costate  equation  Eq.  (11)  has  closed-form  solution 


A(r) 


cos [u;(t  -  r0)]/3 
—  sin[tc;(r  -  t0)]/3 


sin [u(t  -  t0)]/3 
cos[cu(t  -  r0)]/3 


A0  :=  <h(r  -  r0)A0 


where  J3  is  a  3  x  3  unit  matrix.  Define 


(15) 


=  /  lrv(0  cos(uOAT(Qd(  :=  ic(0<K  e  R3 

J  TQ  J  TQ 

Is(t,t0)  =  f  lPv(()sm(u()AT(()d(  :=  f  is(()d(  G  R3 

J  TO  J  TO 


(16) 

(V) 


where  lPv  =  Py/IIPvll-  Note  that  thrust  acceleration  At(-)  is  time- varying,  because  the 
mass  is  changing.  Also  pay  attention  to  the  the  meaning  of  the  time  arguments  in  <h(r  —  r0), 
Jc(-,  •)  and  Js(-,  •)  because  they  will  be  important  in  multi-burn  cases.  Let 


a it  = 


r{r) 

V{r)/u 


Xq  = 


r  o 

Vo/cn 


,  I(t,Tq)  = 


/c(r,  r0) 
Is(t,Tq) 


(18) 


It  can  be  easily  verified  that  the  state  equations  Eqs.  (5)  and  (6)  have  the  following  solution26 

x(t)  =  <L(r  -  r0)a;o  +  T(r)/(r,r0)  (19) 

where 


r(r)  = 


LU 


sin(ct;r)/3  —  cos(car)/3 
cos(car)/3  sin(o;r)J3 


(20) 
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The  thrust  integrals  Ic  and  Is  can  be  evaluated  by  a  quadrature  formula.  For  example,  with 
5  =  (r  —  r0)/ 4,  the  Milne’s  rule  leads  to12 

Ij{Ti  ro)  —  — qq — ~  o)  +  32iJ(r0  +  5)  +  12i,- (r0  +  25) 

+  32iJ(r0  +  35)  +  7ij(ro  +  45)] ,  j  —  c,  s  (21) 

The  values  of  the  primer  vector  pv  at  r0  +  i5  are  needed  in  evaluating  the  thrust  integrals, 
and  they  are  given  by  Eq.  (15)  as  functions  of  A0.  Our  experiences  show  that  with  the  non- 
dimensionalization  described  in  this  report,  the  above  quadratures  are  sufficiently  accurate 
for  powered  flight  up  to  several  hundred  seconds.  Additional  segments  in  time  grid  may  be 
used  for  longer  powered  flight  if  necessary. 

7.5  Solution  for  Coast  Arcs 

For  coast  flight  where  At  =  0,  Ref.26  propagates  the  state  by  continuing  to  use  Eq.  (19) 
which  has  only  the  initial  response  term  now.  We  believe  that  it  is  more  advantageous  to 
use  a  more  accurate  solution  approach  to  the  Keplerian  motion  in  a  coast  arc.  The  /  and 
g  series  in  orbital  mechanics33  are  an  option.  These  series  are  Taylor  series  expansions  of  r 
and  V  in  time,  with  all  the  coefficients  expressed  as  functions  of  r0  and  Vo-  The  inverse- 
square  gravity  Model  is  used,  and  the  accuracy  of  the  solution  is  not  compromised  by  the 
linear  gravity  approximation  even  with  relatively  large  radial  distance  changes.  A  possible 
drawback  of  using  the  /  and  g  series  is  that  for  a  fixed  number  of  terms  in  the  software 
for  the  truncated  Taylor  series,  the  accuracy  deteriorates  as  the  time  of  coast  increases.  Yet 
another  choice  is  to  use  the  well  established  method  by  Goodyear34  to  solve  the  Kepler  initial 
value  problem  as  done  in  Ref.11  This  approach  gives  complete  Keplerian  motion  solution  to 
the  problem  within  the  machine  accuracy  that  includes  the  state  at  a  specified  future  time  as 
the  functions  of  the  current  state  (r0  Vo),  and  the  gradients  of  the  future  state  with  respect 
to  r0  and  Vo ■  This  is  the  method  we  have  adopted  in  this  work. 

The  costate,  on  the  other  hand,  may  continue  to  be  propagated  by  Eq.  (15)  during  the 
coast.  A  second  option  is  to  employ  the  analytical  solution  to  the  costate  in  Keplerian  motion 
developed  by  Lawden.31  This  approach  preserves  the  Model  of  inverse-square  gravity  field, 
but  the  solution  is  given  in  a  rotating  frame  related  to  the  orbital  motion.  Once  the  Kepler 
initial  value  problem  is  solved,  the  corresponding  costate  can  be  readily  transformed  back  to 
the  original  inertial  frame.  In  this  work  we  choose  to  propagate  the  costate  by  Eq.  (15)  for 
its  simplicity.  In  any  approach,  the  state  and  costate  at  the  end  of  the  coast  are  uniquely 
defined  by  their  values  at  the  beginning  of  the  coast  and  the  coast  time. 

7.6  Multiple-Shooting  Formulation 

In  principle,  the  optimal  exo-atmospheric  ascent  problem  from  a  given  initial  condition 
{r0,Vo}  with  a  coast  arc  is  determined  by  8  unknowns:  prg  G  R3,  pVo  G  R3,  toms  (the 
time  when  coast  stops  and  the  second  stage  burn  begins),  and  the  free  final  time  Tj.  The 
corresponding  8  conditions  are  the  /e-terminal  conditions  in  Eq.  (8),  switching  condition  Eq. 
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(14),  transversality  condition  H{jf )  =  0,  and  the  other  6  —  k  transversality  conditions  (see 
Refs.5  and12  for  more  detail  on  how  to  obtain  the  remaining  6  — A;  transversality  conditions). 
Directly  solving  for  the  8  unknowns  to  satisfy  the  8  conditions  is  in  essence  a  single-shooting 
approach.  Difficulties  in  convergence  arise  when  the  coast  time  is  just  Moderately  long 
and/or  the  thrust  acceleration  of  the  last  stage  is  relatively  small,  because  in  these  cases  the 
sensitivity  of  the  problem  increases  and  the  ability  for  the  final  burn  to  steer  the  trajectory 
onto  the  terminal  condition  manifold  weakens. 


To  enhance  the  robustness  of  convergence  of  the  algorithm,  two  additional  nodes  are 
added  to  the  formulation  of  the  numerical  problem.  One  node  is  placed  at  the  end  of  the 
first  powered  stage.  Let  tmeco  >  denote  the  instant  when  the  engine  of  the  first  powered 
stage  shuts  down,  and  recall  that  tmeco  is  considered  specified.  The  other  node  is  naturally 
placed  at  Toms  when  the  coast  ends  and  the  second  powered  stage  begins,  where  tqms  is  to 
be  determined.  Figure  14  gives  an  illustration  of  this  multiple-shooting  formulation.  Define 
the  problem  solution  vector  y  =  col(x  A)  G  R 12 .  In  the  interval  (to,  tmeco),  the  solution 
of  y  is  determined  by  the  condition  y0  at  To  by  Eqs.  (15)  and  (19).  We  will  use  V^jeco 
to  signify  the  value  of  y  at  Tmeco  as  r  — >  Tmeco  from  the  left.  Introduce  two  additional 
to-be-determined  vectors 


Vmeco  ~ 


X  MECO 

\+  >  v oms 

/'MECO 


XOMS 

\  + 

* OMS 


The  propagation  of  y  along  the  coast  arc  in  the  interval  (tmeco,  toms )  is  computed  by  using 
Vmeco  as  ^ie  starting  condition.  In  particular  we  will  denote  such  propagated  value  of  y  at 
toms  by  Voms ■  The  state  and  costate  along  the  last  powered  trajectory  in  ( toms ,  t/)  are 
propagated  from  VqMS ■  The  continuity  of  the  state  and  costate  at  tmeco  and  toms  requires 
the  following  two  conditions  to  be  met 


Vmeco  Vmeco  ~  0 
Voms  ~  Voms  =  0 


(22) 

(23) 


The  above  multiple-shooting  formulation  increases  the  numbers  of  unknowns  by  24  (the 
number  of  scalars  included  in  Vmeco  anfl  Voms)-  The  continuity  conditions  Eqs.  (22)  and 
(23)  provide  the  same  number  of  additional  equations.  Thus  the  dimension  of  the  zero-finding 
problem  is  now  8+24=32.  Note  that  the  solutions  for  the  state  and  costate  make  the  problem 
completely  analytical,  from  function  evaluations  to  Jacobian  computation.  The  evaluation 
of  the  Jacobian  requires  some  care  because  of  the  multiple  segments  of  the  trajectory.  The 
Appendix  provides  several  useful  equations  for  the  part  of  the  Jacobian  involving  the  thrust 
integrals  in  Eqs.  (16)  and  (17).  The  Modest  dimension  of  the  problem  does  not  constitute 
a  heavy  computational  requirement  that  cannot  be  met  more  than  adequately  by  a  desktop 
computer  (see  Section  10  for  some  data  on  CPU  time). 

In  addition  to  the  benefit  of  enhanced  robustness  of  the  algorithm,  another  advantage  of 
this  multiple-shooting  formulation  of  practical  significance  is  the  simplicity  of  the  switching 
condition  in  Eq.  (14)  under  this  setting  (see  Eq.  (29)  later  for  example).  This  feature  further 
contributes  to  improve  the  convergence  on  this  condition.  And  the  complexity  of  the  gradient 
of  the  switching  condition  is  also  reduced  notably,  particularly  so  in  the  case  of  multiple  coast 
arcs,  where  the  switching  condition  will  be  more  complicated.10 
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Figure  14: 


Multiple-shooting  formulation  for  optimal  exo-atmospheric  ascent  with  coast 


7.7  Equality  Constraints 


The  unknowns  in  the  problem  defined  in  Section  6.6  are 

2  —  (Ao,  A MECOi  X~tlECOi  TOMS 1  ^OMS)  XOMS >  Tf )  ^ 


(24) 


In  this  section  we  list  all  the  32  equality  constraints  that  z  must  satisfies  in  order  to  meet 
the  necessary  condition  of  optimal  solution  and  the  continuity  conditions.  The  continuity 


constraints  in  Eqs.  (22)  and  (23)  take  the  explicit  forms  of 

Si{z)  =  §(tmeco  ~  h))Ao  —  A meco  =  0  (25) 

s-2 (z)  =  ${tmeco  ~  to)xq  +  Y{tmeco)I°{jmecoi  to)  —  x^4ECO  =  0  (26) 

s3(z)  =  ®(t0ms  -  tMeco)X\ieco  -  A qms  =  0  (27) 

s4 (z)  =  XCOCLSt{ToMS  ~  TmECOi  ^MBCo)  —  XOMS  =  ^  (28) 


where  the  thrust  integral  I^{tmecOi  ro)  =  (-C  IS)T  has  its  components  defined  as  in  Eqs.  (16- 
17)  with  A0  used  in  the  computation  of  the  primer  vector  pv(-).  The  state  vector  xcoast{roMS~ 
tmecoi  x\/[ECO)  is  the  state  at  Toms  (the  end  of  coast  arc),  propagated  by  the  solution  of 
the  Kepler’s  initial  value  problem  from  tmeco  with  the  initial  condition  x\ieco.  The  scalar 
switching  condition  in  Eq.  (14)  is  conveniently  expressed  as 

s5{z )  —  XqMSXqMS  =  0  (29) 

Let  A f  and  Xf  be  the  final  costate  and  state,  determined  in  [ toms , rf]  through  Eqs.  (15)  and 
(19)  with  the  initial  conditions  and  A  qms  anfl  xqMS,  i.e., 

A  /  =  A(r/,  XqMS)  —  $(t/  —  toms)XqMS  (30) 

xf  =  X(TP  Xoms ,  xoms)  =  $(r/  -  tOMs)x^ms  +  r(r/)/oms+(r/,  tOMS)  (31) 

where  Jom,s+(rj,  toms)  is  the  thrust  integral  defined  in  Eq.  (18)  by  pv(-)  computed  in  the 
interval  [tomSiT/]  with  A  qms  as  the  initial  condition. 
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The  k  orbital  insertion  conditions  plus  the  6 —k  equations  derivable  from  the  transversality 
conditions  (cf.  Refs.11  ,,5  and12)  constitute  another  6  constraints  in  the  form  of 

Sti(z)  —  s6(xf, \f)  —  0  (32) 

The  specifics  of  constraints  in  Eq.  (32)  depends  on  the  particular  set  of  orbital  insertion 
conditions.  References5  and12  provide  some  examples  in  this  regard.  The  last  constraint  is 
the  transversality  condition  on  the  Hamiltonian  H(tj)  =  0  for  the  final  time  Tf  is  free.  This 
condition  can  be  simplified  by  taking  into  consideration  of  the  transversality  condition  on 
Pm(jf)-  Since  m(r/)  is  free,  we  should  have  the  transversality  condition14 

Pm{Tf)  =  0 

Use  of  this  condition  and  the  optimality  condition  of  1  t  =  £V/||fVll  in  H{jf)  =  0  based  on 
Eq.  (10)  leads  to  the  last  constraint 

s7(z)  =  -u2XTfxf  +  At(t f)\\pv(T f)\\  -  -  =  0  (33) 

c 

where  pvijf )  is  the  primer  vector  at  77,  propagated  by  Eq.  (15)  from  toms  with  the  initial 
condition  of  A qms- 


7.8  Numerical  Method  for  Exo- Atmospheric  Ascent 


The  optimal  ascent  problem  in  the  preceding  section  eventually  boils  down  to  a  multivariate 
zero-finding  problem  in  which  a  system  of  nonlinear  algebraic  equations  is  to  be  solved 

s(z)  =0,  z  G  A32  (34) 

where  s(-)  =  (sl5  s2,  s3,  s4,  s5,  s6,  s7)  :  R'i2  — >  R 32  is  the  smooth  vector  function  defined  in 
Section  6.7.  The  classical  Newton- Raphson  method  is  a  common  algorithm  used  for  such 
a  purpose.  While  simple  and  effective  in  many  cases,  the  Newton- Raphson  method  can 
suffer  from  convergence  problems  when  the  initial  guess  is  far  from  the  solution,  and  when 
the  Jacobian  of  the  system  (34)  is  singular  or  nearly  singular.  These  problems  are  handled 
much  better  by  the  renowned  Powell’s  dog-leg  trust-region  method.28  This  method  solves 
the  above  nonlinear  equations  by  minimizing  the  scalar  function  F 

F(z )  =  sT(z)s(z )  (35) 


Obviously  the  solution  to  Eq.  (34)  is  the  solution  to  the  minimum  F .  A  solution  of  local 
minimum  for  F  can  exist  that  is  not  the  solution  of  Eq.  (34)  if  Fmin  ^  0.  But  often  times 
this  is  a  case  where  the  solution  to  Eq.  (34)  does  not  exist  or  the  initial  guess  is  very  poor. 
In  such  a  case  one  would  probably  be  receptive  to  taking  a  solution  that  minimizes  F. 


Powell’s  dog-leg  method  is  in  essence  a  combination  of  the  Newton  and  steepest-descent 
approach.  At  each  iterate  z &  the  update  is  constructed  as 


-i 


s(zk)  +/A 


d  F 
dz 


5z  =  ai 


ds 

dz 
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The  coefficients  ay  and  (3\  are  determined  in  the  algorithm.  When  ay  =  —1  and  =  0,  the 
above  update  reduces  to  the  standard  Newton  update.  This  is  exactly  what  the  algorithm 
does  whenever  the  standard  Newton  update  has  a  size  within  the  trust  region.  Therefore  the 
nice  quadratic  convergence  property  of  the  Newton  method  is  preserved  when  possible.  For 
detail  of  the  algorithm,  the  reader  is  referred  to  Ref.28  In  addition  to  the  well-documented 
superior  performance  and  robustness,  other  attractive  features  of  this  algorithm  include  that 

(i)  exact  Jacobian  of  the  system  Eq.  (34)  is  not  required.  In  fact,  the  Jacobian  can  even  be 
singular  (in  such  a  case  ay  =  0  in  above  update),  a  case  where  the  Newton  method  will  fail; 

(ii)  the  algorithm  will  satisfy  one  of  the  two  stopping  criteria  in  a  finite  number  of  iterations. 
All  these  characteristics  are  appealing  to  our  applications.  We  have  encountered  a  number 
of  cases  with  relatively  challenging  orbital  insertion  conditions  for  which  the  Newton  method 
cannot  converge.  The  code  with  Powell’s  method  on  the  other  hand  finds  the  solution  readily. 

A  singularity  problem  arises  in  the  case  where  the  optimal  coast  time  turns  out  to  be 
zero.  In  such  a  case  the  problem  loses  one  parameter  ( Toms  =  tmeco )  and  the  switching 
condition  in  Eq.  (14),  hence  Eq.  (29),  will  not  be  needed.  With  the  Powell’s  method,  a  simple 
modification  of  the  formulation  can  automatically  adjust  the  problem  to  such  a  case,  so  the 
same  code  can  be  used  for  both  cases  with  and  without  coast  arc.  Define  the  parameterization 

'Tcoastmax  ■  ~  ,  Tcoastmax  "P  2 TmECO  /o r\ 

Toms  = - ^ - sin  2y9  H - - -  (36) 

where  Tcoastmax  is  an  upper  bound  of  the  coast  time.  Clearly  when  P19  =  — vr/2,  the  above 
toms  —  tmeco ,  and  the  coast  time  is  zero.  Now  let  z19  replace  Toms  as  a  parameter  to  be 
found.  The  constraint  in  Eq.  (29)  is  modified  to  be 

«5  =  (1  +  sin  z19)\ omsxoms  =  0  (37) 

When  the  optimal  coast  time=  0,  we  have  sin  iy9  =  —1,  and  S5  =  0  automatically.  The 
problem  reduces  to  a  two-burn  problem,  and  the  Powell’s  method  has  no  difficulty  to  proceed 
to  find  the  solution  without  the  need  to  make  any  changes  to  the  same  burn-coast-burn 
code.  In  contrast,  this  simple  approach  will  cause  difficulty  to  Newton  method,  because 
when  sinS19  =  —  1  (thus  cosiy9  =  0),  the  gradient  of  ds5/dz  is  zero  and  the  Jacobian  of  the 
problem,  ds/dz  e  f?32x32,  is  singular.  As  a  result,  a  separate  code  is  required  to  handle  the 
case  where  the  optimal  coast  is  close  or  equal  to  zero. 
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8 


Analysis  of  the  Condition  H{rt 


=  0 


The  condition  H(jf )  =  0  in  Eq.  (33)  turns  out  to  be  a  very  numerically  difficult  constraint. 
This  difficulty  arises  from  the  fact  that,  unlike  other  constraints  in  Section  6.7,  Eq.  (33) 
can  have  a  great  mismatch  in  numerical  scaling.  Specifically,  the  term  T{jf)/c  is  about  5 
orders  of  magnitude  larger  than  the  rest  of  the  terms  in  Eq.  (33)  when  the  state  and  costate 
vectors  are  dimensionless  and  the  costate  vector  (pj  Py)T  is  chosen  to  have  a  magnitude  on 
the  order  of  unity,  as  usually  done  in  practice.  This  difficulty,  however,  can  be  alleviated  by 
proper  scaling  of  the  problem.  The  governing  equations  for  the  costates  are 


Pr  = 


PV  = 


Pm 


dH 

dr 

dH 

dV 

dH 


dg1 


g0dr 


-Pv 


~Pr 

_ =  T\\pv(t)\\ 

dm  m2(T)g0 


(1) 

(2) 

(3) 


where  the  optimality  condition  1^  =  Iv/lliVll  h&s  been  used  in  Eq.  (3).  Clearly  when  the 
vector  (pj  py  pm)T  is  scaled  by  any  positive  constant,  the  solutions  to  above  costate  equa¬ 
tions  remain  unchanged.  If  the  performance  index  (9)  is  also  scaled  by  the  same  constant,  the 
switching  function  S  in  (12)  will  stay  the  same,  hence  the  optimal  solution  to  the  problem 
will  be  the  same.  Thus,  if  we  choose  to  scale  the  performance  index  and  (pj  py  pm)T  by 
the  value  of  c/T  of  the  second  stage,  and,  with  some  abuse  of  notation  (again),  still  use  pr 
and  pv  to  denote  the  corresponding  vectors  after  the  scaling,  transversality  condition  (33) 
now  becomes 


H(Tf)=p*Vf 


Pvfg(rf)/go  +  AT{Tf)\\pv 


1  =  o 


(4) 


This  version  of  the  condition  H(jf )  =  0  now  admits  costate  vector  ( pJ.(jf )  Py(jf))T  of 
magnitude  comparable  to  those  of  dimensionless  V f  and  r f  (on  the  order  of  unity). 

We  shall  take  the  scaling  analysis  further.  In  fact,  the  above  discussion  also  shows  that  the 
transversality  condition  H{jf)  =  0  is  equivalent  to  the  following  condition  for  any  constant 
k0  >  0 


Prfvf  +  Pvf9(rf)/go  +  AT(Tf)\\pvf  ||  -  k0  =  0 


(5) 


51 


Suppose  for  the  moment  that  the  burn-coast-burn  problem  formulated  in  Section  6  is  solved 
without  explicitly  using  the  condition  H{jf )  =  0.  Rather,  constraint  (33)  is  replaced  by  the 
following  condition 

PTrtV f  +  plfg(r f) / g0  +  AT{Tf)\\pVf\  \  :=  H0(rf)  +  AT{rf)\\pVf\\  >0  (6) 

Once  the  problem  is  solved,  let  k0  =  H0(rf)  +  AT{Tf)\\pVf\\  >  0.  Then,  the  state  and  costate 
solutions  so  obtained  and  the  positive  constant  satisfy  condition  (5).  In  other  words,  the 
costate  vector  (pf  pfl  pm)T  and  the  performance  index  can  be  scaled  by  a  positive  constant 

KqC 

so  that  the  scaled  variables  will  meet  precisely  the  condition  H{jf )  =  0  as  in  the  original 
form  in  Eq.  (33).  This  way,  all  the  necessary  conditions  of  optimality  for  the  problem  are 
now  met.  Let  us  formalize  this  hireling: 

Property  1: 

In  a  free  final-time,  minimum  propellant- consumption,  exo-atmospheric  rocket  flight  prob¬ 
lem,  the  transversality  condition  H{jf )  =  0  is  equivalent  to  the  condition  in  Eq.  (6),  repro¬ 
duced  below  for  the  convenience  of  reference 

Prfvf  +Pvfg(rf)/9o  +  AT(Tf)\\Pvf\\  >  0 

This  property  may  be  of  limited  practical  usefulness  by  itself.  However,  we  shall  show 
next  that  under  two  very  reasonable  assumptions,  condition  (6)  will  be  met  automatically  by 
the  optimal  solution  without  explicit  enforcement.  Consequently,  the  condition  Hfirfl  =  0 
will  be  satisfied  automatically.  The  two  assumptions  that  we  will  use  to  establish  our  result 
are: 

Assumption  1:  The  gravity  held  is  a  Newtonian  inverse-square  force  held. 

Assumption  2:  The  orbital  insertion  conditions  in  Eq.  (8)  are  such  that  they  satisfy  the 
following  condition  in  the  absence  of  thrust  ( T  =  0) 

=  ,(T/)  =  o  (7) 

Assumption  1  means  that  the  unpowered  motion  is  a  Keplerian  orbit.  Assumption  2 
above  is  satished  when  the  constraints  in  0(r/,Vy)  are  expressed  in  terms  of  any  of  the 
Keplerian  orbital  elements  such  as  a  —  aref  =  0  or  e  —  eref  =  0,  where  a  and  e  are  the 
semi-major  axis  and  eccentricity,  respectively,  and  arej  and  eref  are  the  specihed  values  for 
them.  In  fact,  for  such  a  constraint  0'(r)  =  0,  1  <  i  <  k,  for  r  >  ry  because  of  the  constancy 
of  the  orbital  elements.  Other  cases  where  Assumption  2  is  met  include  constraints  on  the 
radius  and  the  magnitude  of  velocity  for  insertion  into  a  circular  orbit,  or  a  non-circular 
(elliptic,  parabolic  or  hyperbolic)  orbit  at  the  periapsis,  or  insertion  into  an  elliptic  orbit 
at  the  apoapsis.  In  the  case  of  non-circular  orbits,  such  a  constraint  meets  the  condition 
( p'flTf )  =  0  only  at  the  orbital  insertion  point  ry.  But  this  is  all  we  need.  In  other  words, 
most  of  typical  orbital  insertion  conditions  can  be  formed  to  satisfy  Assumption  2. 
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With  this  preparation,  we  are  ready  to  formally  make  the  following  claim: 


Property  2: 

Under  Assumptions  1  mid  2,  the  condition 

Ho  (jf)  =  PrfVf  +  Pyfg{rf)/g0  =  0  (8) 

is  always  automatically  satisfied  by  the  solution  to  an  exo- atmospheric  optimal  ascent  prob¬ 
lem. 


Let  us  justify  the  claim.  The  transversality  conditions  on  the  costate  are14 


(  I’r,  \  =  f  OtP'irj,  Vfi/clr,  \ 

V  Pv,  )  V  Vf)/dVf  ) 


(9) 


where  v  G  Rk  is  a  constant  multiplier  vector.  Using  the  expressions  of  pr;  and  pVf  from 
Eq.  (9)  in  iL0(r/),  we  have 


H0(rf)  =  v1 


d4>(rf,Vf )  d  <p(rf,Vf)  .  w 

LVf+  gy  g(rf)/go 


dr  i 


(10) 


Note  that  r' f  —  V f  and  V'f  =  g(rf)/g0  when  T  =  0.  So  in  the  absence  of  thrust,  the  sum 
in  the  above  equation  is  equivalent  to 


Ho(rf)  =  v1 


d<j>(rf,  V f)  ,  d<t>{rf,Vf) 


dr  i 


-r  /  + 


dV, 


V, 


=  <'T<t>\Tf)  =  0 


(11) 


where  Assumption  2  has  been  used  in  arriving  this  condition.  This  conclusion  is  first  arrived 
in  Ref.11  based  on  an  analysis  of  orthogonality  of  the  final  costate  vector  with  respect  to 
the  manifold  of  the  admissible  variations  of  the  final  state  vector.  Since  H{jf )  =  H0(rf)  + 
T{jf)S{Tf )  and  if (t/ )  =  0  is  a  transversality  condition  for  the  optimal  solution,  an  immediate 
corollary  of  =  0  is  that  under  Assumptions  1  and  2,  the  switching  function  satisfies 

S(rf)  =  0  along  the  optimal  solution.  Note  that  S(rf)  =  0  is  not  always  true  because  ry  is 
not  an  interior  switching  point. 

Because  AT{rf)\\pv{Tf)\\  >  0,  condition  (6)  is  always  met  if  ffo(r/)  =  0-  Combining 
Property  1  and  2,  we  see  that  the  enforcement  of  the  transversality  condition  H(tj)  =  0 
is  not  necessary  in  finding  the  solution  to  the  problem  of  present  interest  to  us,  hence  the 
conclusion 


Property  3: 

In  a  free  final-time,  minimum  propellant- consumption,  exo- atmospheric  ascent  problem 
in  which  Assumptions  1  and  2  are  satisfied,  the  transversality  condition  H{jf )  =  0  need  not 
be  enforced  in  the  solution  process. 

Based  on  the  above  conclusion,  one  of  the  to-be-determined  unknowns  in  Eq.  (24)  for 
the  burn-coast-burn  problem  may  be  eliminated  in  correspondence  with  the  removal  of  the 
condition  if(ry)  =  0.  Instead  of  doing  so,  which  has  no  clear  and  convenient  choice,  we  choose 
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to  keep  the  same  unknowns  as  in  Eq.  (24).  But  we  pick  a  simpler  and  easier  replacement  to 
the  condition  H(tj)  =  0.  Numerous  such  choices  exist.  For  instance,  one  could  be 

S7  =  Pv(Tf)PV(Tf )  -  c7  =  o  (12) 

where  c7  >  0  is  a  constant.  The  only  function  of  this  last  constraint  is  to  maintain  the  same 
number  of  equations  as  the  number  of  unknowns  in  Eq.  (24),  and  keep  the  Jacobian  of  the 
constraints  nonsingular.  Therefore,  in  the  numerical  code,  we  do  not  even  need  to  actually 
evaluate  s7  and  attempt  to  enforce  s7  =  0.  Rather,  we  can  always  set  s7  =  0  every  time 
when  the  value  of  s7  is  required  (in  the  case  of  Eq.  (12),  this  would  be  equivalent  to  choosing 
c7  to  be  exactly  equal  to  Pv(rf)PV(Tf)  each  time  s7  is  evaluated).  When  the  linear  gravity 
approximation  is  used,  another  alternative  is  to  use  the  constraint 

S7  =  Pr  ( Tf)pr(Tf )  +  UJ^Py  (t f  )pV  {t f  )  -1=0  (13) 

It  can  be  easily  shown  that  the  magnitude  of  the  vector  A(r)  =  (pj  upy)T  is  constant 
under  the  linear  gravity  Model.  Therefore  if  the  initial  value  of  A(to)  is  scaled  to  have 
unit  length,  the  constraint  in  Eq.  (13)  is  trivially  met.  Such  a  normalizing  constraint  is 
found  to  be  beneficial  to  convergence  when  the  exo-atmospheric  algorithm  interacts  with  the 
endo-atmospheric  algorithm  in  certain  way  as  described  in  Ref.2' 

In  the  end,  we  still  have  the  same  number  of  unknowns  in  Eq.  (24)  as  the  number  of 
constraints  s  =  (si,  S2,  S3,  S4,  S5,  S6,  s7),  where  s7  is  a  trivial  constraint,  and  the  rest  are 
defined  in  the  preceding  section. 

Finally,  as  a  general  interest,  it  is  worthwhile  to  point  out  that  a  similar  claim  can  be 
made  for  the  minimum-time  optimal  ascent  problem  for  insertion  into  a  Keplerian  orbit.  In 
such  a  problem,  the  transversality  condition  on  the  Hamiltonian  is14 

h(t,)  =  pTr,V,  +  PTVlg(r,)/g„  +  =  1 

Suppose  that  the  problem  is  solved  without  explicitly  enforcing  the  above  condition.  Prop¬ 
erty  2  above  still  holds,  which  together  with  1  t(t/)  =  Pv(rf) /\\Pv(rf)\\  leads  to 

H(Tt)  =  n  T/)IIPt(?)l1  >0  (14) 

m(Tf)g0 

We  now  make  note  of  two  observations  of  invariance  pertinent  to  the  problem: 

1.  The  costate  vector  (pj  Py)T  can  be  scaled  by  any  constant  and  the  costate  equa¬ 
tions  (1-2)  are  still  satisfied; 

2.  When  the  costate  vector  (pj  Py)T  is  scaled  by  any  positive  constant,  the  resultant 
optimal  control  from  the  optimality  condition  and  the  state  trajectory  will  remain  the 
same. 


Hence  if  we  scale  the  costate  by  a  positive  constant 


«2 


™(.Tf)9o 

T(rf)\\Pvf\\ 
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the  optimal  solution  will  remain  unchanged,  and  we  will  now  have  H(tj)  =  1  without 
explicitly  enforcing  it.  Indeed,  the  maximum  principle  by  Pontryagin  et  a/35  only  requires 
H{jf )  >  0  for  the  minimum-time  problem,  which  is  already  met  by  Eq.  (14). 
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9 


Orbital  Insertion  Modes 


As  noted  a  few  times  earlier  in  the  report,  the  orbital  insertion  conditions  are  specified  in 
terms  of  k  equality  constraints  on  the  final  state;  More  specifically,  typical  inertial  insertion 
conditions  are  defined  by  a  subset  of  the  six  target  orbital  elements.  This  section  will  describe 
several  terminal  Modes  that  specify  various  commonly  seen  launch  missions.  Another  com¬ 
ponent  of  this  section  is  to  show  how  the  Lagrange  multipliers  in  the  transversality  conditions 
of  the  costate  may  be  eliminated  by  manipulating  the  transversality  conditions  and  orbital 
insertion  conditions  cleverly.  The  elimination  of  these  unknown  multipliers  has  two  major 
advantages:  first  it  reduces  the  number  of  unknowns  needed  to  be  found;  more  importantly, 
it  increases  the  algorithm  convergence  reliability  significantly.  The  second  point  arises  from 
the  fact  that  these  multipliers  have  no  physical  meaning,  and  they  may  differ  in  magnitude 
with  that  of  the  rest  of  variables  considerably.  Having  to  solve  for  these  multipliers  would 
make  the  solution  process  more  difficult  in  general. 

In  this  section  the  final  state  and  costate  are  denoted  by 

*' =  (  V,  )  (1) 

*=(£;)  (2) 


9.1  Mode  31 

This  Mode  consists  of  three  orbital  insertion  conditions  defined  by  the  desired  target  orbital 
element  values  for  the  semi- major  axis  a*,  eccentricity  e* ,  and  inclination  i*.  The  ascending 
node,  true  anomaly,  and  argument  periapsis  are  considered  free.  Since  the  point  of  insertion 
into  the  target  orbit  is  not  specified,  this  is  a  so-called  free  attachment  Mode,  and  the  final 
flight  path  angle  7/  is  unconstrained.  When  the  insertion  flight  path  angle  is  specified  such 
as  7 f  =  0,  thereby  requiring  the  insertion  point  be  at  either  the  perigee  or  apogee  of  an 
elliptic  orbit,  one  additional  constraint  will  be  needed.  (See  Mode  43/44  below). 

It  should  be  noted  that  naming  convention  for  the  orbital  insertion  Modes  is  such  that  the 
first  digit  (“3”  in  this  case)  indicates  the  number  of  orbital  insertion  constraints  in  each  case. 
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The  second  digit  simply  designate  different  different  variation  of  orbital  insertion  conditions 
while  the  number  of  constraints  is  the  same. 


A  unit  vector  is  constructed  that  relates  the  vehicle’s  final  position  and  velocity  to  the 
target  orbit  inclination.  Define  a  z-axis  unit  vector  in  the  Earth  Centered  Inertial  (ECI) 
frame  as 

"  0 
i 

ljV£  = 


0 
1 

This  vector  transformed  to  the  Guidance  frame  is  given  by 

1 NG  =  Tep  InE 

where  Tep  is  the  transformation  matrix  as  defined  in  section  2.4. 


(3) 


(4) 


From  orbital  mechanics,  the  magnitude  of  the  required  final  angular  momentum  vector 
of  the  vehicle  at  the  orbital  insertion  point  is  given  by 

h*  =  \J  a*(l  —  e *2)  (5) 

Using  the  above  expressions,  the  three  orbital  insertion  conditions  T'  G  R3  can  be  written 


as 


5  (<•/  ^  V, )T  (r /  x  V,)  -  i ft'2  =  0 

1  1  1 

-VfTVf  -  —  +  —  =  0 

2  1  rf  2a* 


1  ngT  i?f  x  V f)  —  h*  cos  i*  —  0 

where  ry  =  ||r^||. 

Additional  Terminal  Conditions  After  Elimination  of  Lagrange  Multipliers 


(6) 

(7) 

(8) 


The  6  transversality  equations  in  Eqs.  (25-26)  in  this  contain  3  Lagrange  multipliers 

{v  e  R3): 

1 

Pr;  ~  [V f  x  hf)  vx  -  — r/i/ 2  -  (Vf  x  1NG)  u3  =  0  (9) 

PVf  +  (rf  x  hf)  -  V fv2  +  {rfx  1NG)  u3  =  0  (10) 

After  using  the  conditions  (6-8)  and  the  transversality  equations  to  eliminate  the  v  (the 
algebraic  manipulations  are  lengthy  and  thus  not  shown  here),  three  additional  terminal 
conditions  are  resulted  in: 

(■ hfTPrf )  [ hfT  {rf  x  Ijvg)]  +  ( hfTPvf )  [ hfT(vf  x  !jvg)]  =  0  (11) 

rf3VfTPrf-rfTpV/=  0  (12) 
PrfT  irf  xVN)  (r /  xhf)T(Vfx  rN)  +  pVfT  ( V fxrN)(Vfx  hf)T  (rf  x  V N)  =  0  (13) 
where  hf  =  rf  xV  f  is  the  angular  momentum  vector,  rN=VfXlNG,  and  V  N  =  V  fXlNG. 
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9.2  Mode  41 


Mode  41  consists  of  four  orbital  insertion  conditions  defined  by  the  desired  target  orbital 
element  values  for  the  semi-major  axis  a*,  eccentricity  e*,  inclination  i*,  and  ascending 
node  fl*.  The  true  anomaly  and  argument  periapsis  are  considered  free.  This  is  also  a  free 
attachment  Mode  in  that  the  final  flight  path  angle  7/  is  unconstrained  allowing  the  insertion 
point  to  be  located  at  any  point  on  the  target  orbit. 

A  unit  angular  momentum  vector  of  the  desired  target  orbital  plane  is  constructed  from 
the  target  inclination  and  ascending  node.  This  vector,  defined  in  the  ECI  frame  is  given  as 


1  he  — 


sin  Q*  sin  i* 
—  cos  Q*  sin  i* 
cos  i* 


(14) 


This  vector  transformed  to  the  Guidance  frame  is  given  by 


1 HG  —  Tep  1  HE  (15) 

Again,  the  magnitude  of  the  required  final  angular  momentum  vector  of  the  vehicle  at 
the  orbital  insertion  point  is  given  by 

h*  =  \J  a*(  1  —  e*2)  (16) 

Using  the  above  expressions,  the  four  orbital  insertion  conditions  $  6  U4  can  be  written 
as 


(17) 

(18) 

where  ry  =  ||ry||  and  hf*  —  h*luG 

Additional  Terminal  Conditions  After  Elimination  of  Lagrange  Multipliers 


O  VfTVf 


rfxVf- 
1 

rf  +  2  a 


hf*  —  0 

1 


g  ir 


=  o 


The  6  transversality  conditions  involving  the  Lagrange  multiplier  vector  u  G  R4  are  now 
given  by 


prf  +  skew  (V f) 


v\ 

u2 

"3 


~  yi’’/1'!  =  0 


pv  -  skew  (77) 


v\ 

V3 


—  V  fV4  =  0 


(19) 

(20) 


where  skew( V  f )  and  skew(r f)  are  the  3x3  skew  matrices  formed  by  the  elements  of  V f  and 
rf,  respectively.  After  much  algebraic  operation  to  eliminate  u  from  the  above  equations, 
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we  can  have  two  additional  terminal  conditions  free  of  v. 


Pr/Vf  -  7fPv,Trf  =  0 

(21) 

PrfT(rf  x  hf*)+pVfT(Vf  x  hf*)  =  0 

(22) 

9.3  Mode  43/44 

This  Mode  involves  four  orbital  insertion  constraints  defined  by  the  desired  target  orbital 
element  values  for  the  semi- major  axis  a*,  eccentricity  e*  inclination,  i*7  and  a  final  flight 
path  angle  7/,  which  is  set  to  zero.  Therefore  the  insertion  point  is  either  at  the  perigee 
(Mode  43)  or  apogee  of  the  target  orbit  (Mode  44)  and  the  final  dimensionless  position  and 
velocity  are  given  by 

(23) 

(24) 


r}=a"(  1  —  e*) 


Vf  = 


1  +  e* 


for  the  perigee  and 


a*(l  —  e* 
r}  =  a*(  1  +  e*) 


Vf  = 


a*(l  +  e*) 


(25) 

(26) 


for  the  apogee. 


A  unit  vector  is  constructed  that  relates  the  vehicle’s  final  position  and  velocity  to  the 
target  orbit  inclination.  Define  a  z-axis  unit  vector  in  the  ECI  frame  as 


Iace  — 


0 

0 

1 


(27) 


This  vector  transformed  to  the  Guidance  frame  is  given  by 

IjVG  =  Tep  Iwe 


(28) 


Using  the  above  expressions,  the  four  orbital  insertion  conditions  G  R4  can  be  written 
as 


-VfTVf  -  -Vf  =  0 

2  /  1  2  f 


1  ngt  (rf  xVf)  —  r'fVj  cos**  —  0 

rfTVf=0 


(29) 

(30) 

(31) 

(32) 


59 


Additional  Terminal  Conditions  After  Elimination  of  Lagrange  Multipliers 


The  6  transversality  conditions  involving  the  Lagrange  multiplier  vector  v  e  R4  are  now 
given  by 


Prf  -rfu I  -  (Vf  x  1NG)  v3-Vfv 4  =  0  (33) 

PVf  ~V fu2+  (rf  X  lNG)u3  -  rfU4  =  0  (34) 

Performing  strategic  dot  and  cross  product  operations  on  these  conditions  enables  the 
elimination  of  the  multiplier  vector  v  e  R4.  This  results  in  2  additional  independent  terminal 
constraint  equations  free  of  the  multipliers  given  by 

VJTPrfrf-rfTpVfV?  =  0  (35) 


(rf  x  Vf)Tprf 
+ 


( rf  X  Vf)T(rf  X  Ijvg) 

( rf  x  V f)T pVf 

(■ rf  X 

V,f  ( V ,  x  1NG)' 

0 


(36) 


9.4  Mode  46 


This  Mode  involves  four  insertion  conditions  defined  by  the  desired  final  position  magnitude 
rjb  velocity  magnitude  Vf,  inclination  i* ,  and  flight  path  angle  7^  which  is  not  necessarily 
zero.  If  values  of  r*j,  Vf,  and  7))  are  chosen  properly,  the  orbital  insertion  point  can  be  placed 
anywhere  in  the  orbit  using  this  Mode.  So  Mode  43/44  is  a  special  case  to  Mode  46. 

Again,  a  unit  vector  is  constructed  that  relates  the  vehicle’s  final  position  and  velocity 
to  the  desired  inclination.  Define  a  z-axis  unit  vector  in  the  ECI  frame  as 


ljV.E  — 


0 

0 

1 


This  vector  transformed  to  the  Guidance  frame  is  given  by 

liVG  =  Tep  Ine 


(37) 


(38) 


The  four  resulting  insertion  conditions  can  be  written  as 

It  1  *2 

2r!  rf  ~  2 rl  =  0 

\viTVf  -  \vp  =  0 

1  ngT  {f'f  x  V f)  —  r*fVf  cos  7))  cos  i*  =  0 
r/1  V f  —  r*fVf  sin  7/  =  0 


(39) 

(40) 

(41) 

(42) 
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Additional  Terminal  Conditions  After  Elimination  of  Lagrange  Multipliers 


The  6  transversality  cantons  involving  the  Lagrange  multiplier  vector  u  E  R4  are  now 
given  by 


Prf  ~  rfv i  -  (V f  x  1NG)  v3  -  V fuA  =  0 


(43) 

(44) 


PVf  ~  v fu 2  +  (rf  x  Ijvg)  u3  -  rfu4  =  0 
The  2  final  manipulated  additional  terminal  conditions  free  of  the  Lagrange  multipliers  are 


y f  Prf)  r}  -  [rf  pVf  )  V )*  +  r*fv;  sin 7}  Vf  pVf  -  rf  p 


rf 


=  0 


(45) 

/>■/>,,)  (hfTrN)  +  (h/Pv,)  (h/TVN)  =  0  (46) 

where  hf  —  r  fXV  f  is  the  angular  momentum  vector,  r  N  =  rj  x  1  no,  and  V  N  =  V  f  xI^g- 

9.5  Mode  51 

Mode  51  consists  of  five  orbital  insertion  conditions  defined  by  the  desired  target  orbital 
element  values  for  the  semi-major  axis  a*,  eccentricity  e*,  inclination  i*,  ascending  node  0*, 
and  and  a  final  flight  path  angle  7/  =  0.  This  is  the  Mode  to  be  used  for  a  launch-to- 
rendezvous  mission.  The  insertion  point  is  at  the  perigee  or  apogee  of  the  target  orbit  and 
the  final  position  and  velocity  are  given  again  as 


r}  =  a*(l-e*) 


V*  — 
vf  ~ 


1  +  e* 


for  the  perigee  and 


a*(l  —  e* 
r*f  =  a*(l  +  e*) 


Vf  = 


a*(l  +  e*) 


(47) 

(48) 

(49) 

(50) 


for  the  apogee.  As  was  done  for  Mode  41,  a  unit  angular  momentum  vector  of  the  desired 
target  orbital  plane  is  constructed  from  the  target  inclination  and  ascending  node  and  is 
given  as 

sin  Q*  sin  i* 

1  he  —  I  —  cos  Q*  sin  i*  (51) 

cosi* 


This  vector  transformed  to  the  Guidance  frame  is  given  by 

1  hg  —  Tep  1-he 


(52) 
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Again,  the  magnitude  of  the  required  final  angular  momentum  vector  of  the  vehicle  at 
the  orbital  insertion  point  is  given  by 


h*  =  \J a*(l  -  e*2)  (53) 

Using  the  above  expressions,  the  hve  orbital  insertion  conditions  \I/  G  R5  can  be  written 
as 


■Tr,~\rf  =  0 

(54) 

v,  -  ip*2  =  0 

(55) 

rfTVf  =  0 

(56) 

rfT  Ihg  =  0 

(57) 

V /T1hg  =  0 

(58) 

Additional  Terminal  Conditions  After  Elimination  of  Lagrange  Multipliers 

The  6  transversality  conditions  involving  the  Lagrange  multiplier  vector  v  E  R5  are  now 
given  by 


Prf  ~  r fvi  —  V fu3  -  lHGu4  =  0  (59) 

Pvf  ~  V f"2  -  rfis3  -  1Hgv5  =  0  (60) 

The  final  manipulated  additional  terminal  condition  free  of  the  Lagrange  multipliers  is 
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10 


Combining  Atmospheric  and  Vacuum 
Algorithms 


The  endo-atmospheric  and  exo- atmospheric  ascent  guidance  algorithms  developed  in  the 
preceding  parts  of  this  report  must  work  together  to  produce  a  complete  optimal  ascent 
trajectory  from  liftoff  to  orbital  insertion.  The  challenge  in  this  approach  of  combining 
the  atmospheric  ascent  algorithm  with  the  vacuum  ascent  algorithm  lies  in  the  integration, 
since  the  vacuum  algorithm  is  not  just  a  special  case  of  the  atmospheric  algorithm  with  zero 
atmospheric  density.  Instead  of  solving  the  entire  problem  with  one  algorithm,  each  of  the 
algorithms  solves  one  part  of  an  optimal  control  problem  and  the  trajectories  in  each  part  are 
then  pieced  together.  This  division,  however,  poses  a  unique  problem  for  the  atmospheric 
ascent  algorithm  because  the  end  conditions  for  the  atmospheric  portion  of  the  ascent  are 
not  defined  a  priori.  Rather,  they  are  a  result  of  the  final  solution.  For  instance,  the  usual 
multiple-shooting  formulation  at  the  junction  point  will  have  difficulty  since  two  different 
algorithms  are  at  work,  each  covering  only  one  part  of  the  entire  trajectory. 

Our  approach  is  to  integrate  the  algorithms  through  the  iteration  on  the  state  at  the  end 
point  Ti  of  the  first  stage  where  the  atmospheric  algorithm  stops  and  the  vacuum  algorithm 
begins.  Figure  15  illustrates  the  point.  From  a  starting  estimate  of  the  state  at  Ti,  the  vacuum 
algorithm  generates  an  optimal  trajectory  to  the  target  orbit  and  returns  the  corresponding 
costate  at  T\.  The  atmospheric  algorithm  in  turn  solves  a  TPBVP  in  which  the  initial  state 
is  that  at  the  lift-off  and  the  final  costate  is  the  one  just  found  by  the  vacuum  algorithm  at 
T\.  This  atmospheric  ascent  solution  will  provide  a  new  state  at  T\.  The  above  process  then 
repeats  until  the  states  found  at  Ti  in  two  consecutive  cycles  are  practically  the  same. 

The  algorithmic  implementation  of  the  above  approach  is  given  below.  Starting  from  an 
initial  guess  of  the  state  at  Ti,  denoted  by  ,  our  integration  algorithm  proceeds  with  the 
iterations  as  follows: 

1.  Set  k  =  0. 

2.  With  the  known  x^\  the  vacuum  algorithm  finds  the  optimal  ascent  trajectory  in 
(ti,  Tf ),  possibly  including  the  coast,  from  x ^  to  the  targeting  condition  Eq.  (8).  As 
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Figure  15:  Integration  of  endo-  and  exo-atmospheric  optimal  trajectories 


(k) 

a  part  of  the  solution  the  vacuum  algorithm  returns  the  corresponding  costate  pT1  at 
n- 

3.  The  just  obtained  p^  is  used  as  the  required  boundary  condition  at  T\  for  the  costate 
in  the  interval  (0,  T\)  for  the  atmospheric  algorithm.  The  initial  state  at  r  =  0  is 
known.  The  TPBVP  for  the  atmospheric  algorithm  is  then  well  defined.  Upon  the 
convergence  of  the  algorithm,  a  new  state  at  T\  is  found  as  a  result,  and  this  new  state 
is  denoted  as  cc^+P. 

4.  If  ||£C^+P  —  ami'll  <  5  for  some  pre-selected  small  constant  5  >  0,  set  xT] 
stop.  Otherwise,  let 

x%+1)  =  ex^+1)  +  (1  -  e)®^ 

where  0  <  £  <  1  is  a  constant.  Set  k  =  k  +  1  and  return  to  Step  2  above 

At  the  conclusion  of  the  above  process,  both  the  state  x  and  costate  p  are  continuous  at 
T\ .  All  the  necessary  conditions  for  the  optimal  ascent  problem  are  satisfied  in  [0,  r/].  In 
other  words,  the  solutions  obtained  for  :»(•)  and  p(-)  in  the  two  algorithms  form  a  continuous 
extremal  for  the  complete  optimal  ascent  problem. 

The  above  search  for  the  correct  state  xT1  essentially  constitutes  a  fixed-point  problem. 
Consider  the  case  when  £  =  1  in  Eq.  (1).  We  may  represent  the  process  of  obtaining  ai^+P 

based  on  Xt\}  by  the  mapping  x^+l^  =  M.{x^).  The  above  procedure  therefore  amounts  to 
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finding  the  solution  to  the  fixed-point  equation 

xT1  =  M(xn)  (2) 

In  this  context,  Eq.  (1)  with  £  ^  1  simply  is  a  modified  fixed-point  iteration  in  the  form  of 

x(k+1)  =  eM(xW)  +  (1  -  e)® W  (3) 

Evidently  the  fixed-point  solution  to  Eq.  (3)  for  any  £  ^  0  is  the  solution  to  Eq.  (2)  as  well. 
The  introduction  of  an  appropriate  £  in  Eq.  (3)  can  aid  the  convergence  of  the  fixed-point 
iteration  which  may  not  always  happen  with  £  =  1.  A  similar  technique  is  discussed  in  Ref.37 
In  our  application  this  technique  has  a  geometric  interpretation:  an  £  <  1  in  the  early  cycles 
of  the  iterations  can  prevent  xl^1^  from  departing  too  far  from  x^\  an  concurrence  which 
in  some  cases  may  result  in  the  trajectory  going  through  the  Earth.  Once  a^+1)  and  Xr\} 
are  sufficiently  close,  an  e  =  1  can  be  used  to  further  speed  up  the  convergence. 

In  fact,  the  need  to  solve  the  special  endo- atmospheric  TPBVP  in  Step  2  above  inspires 
another  fixed-point  formulation  of  this  sub-problem,  and  consequently  a  fixed-point  iteration 
algorithm  may  be  used  to  solve  this  particular  problem.  Appendix  III  provides  some  further 
exposition  on  this  aspect. 
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V&V  of  Exo- Atmospheric  Ascent 
Guidance  Algorithm 


The  objective  in  this  section  is  to  verify  and  validate  the  analytical  multiple  shooting  (AMS) 
exo- atmospheric  ascent  guidance  algorithm  developed  here  by  comparing  the  solutions  from 
the  algorithm  to  those  obtained  by  using  an  industry  standard  trajectory  optimization  soft¬ 
ware.  Toward  this  purpose,  41  different  mission  scenarios  are  designed  and  tested  with  a 
two-stage  launch  vehicle  Model.  All  the  orbital  insertion  Modes  in  Section  9  are  tested.  To 
verify  the  results  found  by  the  AMS  algorithm,  an  industry  standard  aerospace  trajectory 
optimization  software,  called  Optimal  Trajectories  by  Implicit  Simulation29  (OTIS),  is  em¬ 
ployed  to  compute  the  burn-coast-burn  trajectory  under  the  identical  condition.  Further 
comparison  is  done  with  closed-loop  simulated  trajectories.  In  the  closed-loop  simulations, 
the  AMS  algorithm  serves  as  the  guidance  algorithm,  generating  the  optimal  exo-atmospheric 
trajectory  from  the  current  state  to  the  orbital  insertion  point  at  a  guidance  cycle  of  1  Hz. 
The  actual  trajectory  is  simulated  by  integrating  Eqs.  (1-3)  with  an  inverse-square  gravity 
field.  The  optimal  thrust  direction  and  throttle  command  (full  thrust  or  coast)  are  deter¬ 
mined  by  the  AMS  algorithm.  The  closed-loop  simulations  are  an  ultimate  check  for  the 
validity  of  the  open-loop  solution  found.  If  the  closed-loop  trajectory  matches  the  open-loop 
one  closely,  the  approximations  adopted  in  obtaining  the  open-loop  solution  are  well  justified. 


11.1  Two-Stage  Launch  Vehicle 

A  two-stage  vehicle  configuration  is  used  for  the  verification  and  validation  of  the  analytical 
multiple  shooting  (AMS)  algorithm  for  burn- coast-burn  vacuum  ascent  trajectory  optimiza¬ 
tion  problem.  The  vehicle’s  first  stage  is  the  so-called  “Super  X-33” .  This  is  a  vehicle  Model 
with  all  the  data  identical  to  those  of  the  X-33,  except  that  the  specific  impulse  of  the  engine 
is  doubled.  The  second  stage  is  the  X-37  piggy-backed  on  the  Super  X-33.  The  launch 
site  for  all  test  cases  shown  below  is  at  Kennedy  Space  Center  (KSC),  (longitude=— 80.85°, 
latitude=28.29°).  This  information  is  necessary  in  establishing  the  Guidance  inertial  frame 
for  the  launch  as  well  as  the  required  Guidance  to  ECI  coordinate  transformation.  The  com¬ 
putation  of  the  optimal  first  stage  ascent  trajectory  through  the  atmosphere  is  divided  into 
two  segments:  the  endo-atmospheric  portion  and  the  exo-atmospheric  portion,  separated  at 
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a  given  altitude  around  90  km.  The  exo-atmospheric  portion  of  the  Super  X-33  is  regarded 
as  the  first  upper  (vacuum)  stage.  After  the  burn-out  of  the  X-33,  the  X-37  continues  the 
optimal  ascent  trajectory  by  first  coasting  for  an  optimal  duration  and  then  starting  burn 
of  the  X-37  until  orbital  insertion.  Data  used  in  the  AMS  algorithm  for  the  Super  X-33  and 
X-37  are  given  below. 


Table  4:  Vehicle  data  of  the  first  vacuum)  stage  of  “Super  X-33” 


Thrust  (vac) 

2,303,000  N 

hp  (vac) 

711.395  sec 

Empty  weight 

37,557  kg 

Propellant  weight 

27,000  kg 

Table  5:  X-37  Data 

Thrust  (vac) 

29,269  N 

hP  (vac) 

330  sec 

Empty  Weight 

1,270  kg 

Propellant  Weight 

3,628.74  kg 

Payload  Weight 

544.31  kg 

11.2  Mode  31  /  43  Comparisons 

As  mentioned  previously,  Mode  31  consists  of  three  orbital  insertion  conditions  defined  by 
the  desired  target  orbital  element  values  for  the  semi-major  axis  a*,  eccentricity  e* ,  and 
inclination  i*.  This  is  a  free  attachment  Mode  in  that  the  final  flight  path  angle  7/  is  not 
constrained.  Mode  43  requires  the  same  orbital  insertion  conditions  as  Mode  31  in  addition 
to  constraining  the  final  flight  path  angle  to  zero  such  that  the  insertion  point  is  the  perigee 
of  the  target  orbit.  A  set  of  test  cases  with  target  orbital  insertion  conditions  as  well  as  an 
identical  set  for  Mode  43  including  the  final  7/  =  0  constraint  is  given  below  in  Table  6  .  In 
addition,  three  circular  target  orbits  for  Mode  43  arc  listed. 

Results  for  Mode  31  and  43  test  cases  in  Table  6  are  listed  in  Table  7  and  8  respec¬ 
tively.  Recalling  that  the  burn  time  of  the  first  vacuum  burn  is  fixed  and  determined  from 
propellant  availability,  the  difference  in  converged  second  burn  time  directly  translates  into 
the  difference  in  deliverable  payload  mass,  or  vehicle  performance.  Both  the  second  vacuum 
burn  time  and  mass  are  listed  for  convenience.  In  terms  of  vehicle  performance,  all  methods 
match  quite  closely  and  is  very  noteworthy  in  that  OTIS  is  a  direct  approach  using  colloca¬ 
tion,  a  very  different  method  than  the  indirect  AMS  method  that  finds  the  trajectory  based 
on  optimal  control  necessary  conditions.  Further,  the  extent  to  which  the  open-loop  AMS 
trajectory  matches  the  full  non-linear  gravity  closed-loop  simulation  driven  by  a  current  con¬ 
dition  open-loop  AMS  solution  at  each  guidance  cycle  is  reassuring.  The  small  differences 
in  the  open-loop  AMS  solution  and  the  converged  closed-loop  simulation  is  to  be  expected 
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Tabic  6:  Mode  31  fe  43  Target  Orbits 


Case 

Mode 

perigee  alt.  (km) 

e 

i  (deg) 

7/  (deg) 

1 

31 

300 

0.1 

51.6 

Free 

2 

31 

500 

0.1 

51.6 

Free 

3 

31 

1000 

0.1 

51.6 

Free 

4 

31 

300 

0.3 

51.6 

Free 

5 

31 

500 

0.3 

51.6 

Free 

6 

31 

1000 

0.3 

51.6 

Free 

7 

31 

500 

0.1 

28.5 

Free 

8 

31 

500 

0.3 

28.5 

Free 

9 

43 

300 

0.1 

51.6 

0.0 

10 

43 

500 

0.1 

51.6 

0.0 

11 

43 

1000 

0.1 

51.6 

0.0 

12 

43 

300 

0.3 

51.6 

0.0 

13 

43 

500 

0.3 

51.6 

0.0 

14 

43 

1000 

0.3 

51.6 

0.0 

15 

43 

500 

0.1 

28.5 

0.0 

16 

43 

500 

0.3 

28.5 

0.0 

17 

43 

1000 

0.0 

51.6 

0.0 

18 

43 

500 

0.0 

51.6 

0.0 

19 

43 

500 

0.0 

51.6 

0.0 

and  results  from  approximations  made  in  the  Analytical  Multiple-Shooting  problem  formu¬ 
lation  such  as  the  linear  gravity  approximation  and  thrust  quadrature  approximations  for 
each  burn.  During  the  coast  phase  the  launch  vehicle  is  trading  kinetic  energy  for  potential 
energy.  The  higher  the  orbital  insertion  altitude,  the  longer  the  coast  arc  tends  to  be,  and 
the  more  critical  is  its  length  to  the  performance  of  the  vehicle.  Figures  16  and  17  show  the 
closed-loop  altitude  and  velocity  profiles  along  with  their  corresponding  body  axis  pitch  and 
yaw  angle  time  histories  with  respect  to  the  guidance  frame  defined  with  a  2-3-1  rotation 
sequence.  The  small  yaw  angles  in  figure  17  are  due  to  the  fact  that  the  launch  azimuth  is 
oriented  in  the  direction  defined  by  the  target  orbit  inclination  so  as  to  minimize  the  yaw 
maneuvers,  a  common  practice  in  launch  vehicle  guidance.  Some  cases  were  found  where, 
depending  on  insertion  conditions,  higher  altitudes  did  not  always  translate  to  a  longer  coast 
arc,  however,  the  presence  of  the  optimized  coast  arc  remains  very  significant  to  vehicle  per¬ 
formance. 

The  performance  index  on  the  final  mass  in  the  optimal  burn-coast-burn  problem  appears 
to  be  fairly  insensitive  to  the  coast  time.  Indeed,  significant  differences  in  the  converged  coast 
time  between  OTIS  and  the  AMS  method  does  not  result  in  significant  differences  in  orbital 
insertion  mass.  This  disagreement  in  optimal  coast  times  becomes  much  more  noticeable  for 
higher  altitude  orbital  insertions.  The  two  methods  do,  however,  find  similar  solutions  for 
the  lower  altitude  orbits.  These  results  are  illustrated  in  Figures  18  and  19,  which  show  the 
altitude  and  velocity  profiles  for  both  OTIS  and  the  closed-loop  AMS  simulation  for  cases 
2  and  3.  The  small  changes  in  insertion  mass  even  for  disproportionately  large  changes  in 


coast  time  results  in  a  “flat  optimum”  and  it  was  found  that  direct  methods  such  as  OTIS 
tend  to  converge  to  different  solutions,  most  notably  in  the  converged  coast  time,  depending 
on  initial  guess  inputs.  The  AMS  method,  an  indirect  method,  did  not  suffer  from  this 
problem  and  was  consistently  able  to  converge  to  the  same  solution  for  the  same  case  for 
large  changes  in  the  initial  guess. 

An  important  comparison  between  all  Mode  31  and  the  respective  Mode  43  test  cases  is 
seen  from  the  value  of  the  true  anomaly  at  insertion.  Mode  43  requires  the  insertion  point 
be  the  perigee  of  the  resulting  target  orbit,  which  evident  from  the  table  is  not  always  the 
optimal  insertion  point.  This  result  is  consistent  with  both  the  AMS  method  and  OTIS. 
This  is  exemplified  by  test  case  6  and  14  which  results  in  a  difference  in  second  burn  time  of 
approximately  5  seconds  equating  to  a  deliverable  payload  difference  of  45.25  kg,  a  substantial 
difference  assuming  payload  mass  of  500  kg.  Figures  20  and  21  show  the  optimized  ascent 
trajectory  and  insertion  target  orbit  for  case  6.  This  illustrates  the  non-perigee  optimal 
insertion  point  and  shows  the  true  anomaly  at  insertion.  Figure  21  is  viewed  in  the  direction 
of  the  target  orbit  angular  momentum  vector  allowing  this  to  be  seen  easily. 


Table  7:  Mode  31  Results 


Case 

Method 

Coast  (sec) 

2nd  Burn  (sec) 

Final  Mass  (kg) 

True  Anomaly  (deg) 

OTIS 

24.6240 

334.0512 

2421.7212 

7.43432 

1 

AMS 

20.0371 

334.3829 

2416.8352 

0.00000 

AMS  closed-loop 

26.4927 

331.7147 

2438.4000 

6.21815 

OTIS 

135.2167 

345.9007 

2314.5505 

11.41668 

2 

AMS 

126.1379 

346.8496 

2304.0113 

0.00000 

AMS  closed-loop 

134.5228 

343.6698 

2329.8000 

11.13721 

OTIS 

293.8752 

389.6641 

1918.7385 

10.29114 

3 

AMS 

226.0024 

390.5117 

1908.8694 

0.00000 

AMS  closed-loop 

243.9442 

388.0392 

1922.5500 

12.08240 

OTIS 

54.2115 

386.7927 

1944.7088 

2.38447 

4 

AMS 

43.9989 

386.6586 

1943.7401 

0.00000 

AMS  closed-loop 

53.1648 

384.5900 

1958.7500 

1.29348 

OTIS 

152.4945 

395.3354 

1867.4451 

6.08704 

5 

AMS 

132.0432 

396.0864 

1858.4179 

0.00000 

AMS  closed-loop 

142.8387 

393.3418 

1877.3000 

6.75834 

OTIS 

257.3082 

427.3834 

1577.5918 

11.21645 

6 

AMS 

207.8629 

429.6897 

1554.3082 

0.00000 

AMS  closed-loop 

226.6251 

425.4743 

1587.7000 

12.72862 

OTIS 

129.0704 

345.5416 

2317.7982 

12.07388 

7 

AMS 

129.1712 

345.9306 

2312.3279 

0.00000 

AMS  closed-loop 

137.6387 

342.7682 

2338.8500 

10.95298 

OTIS 

126.4071 

394.4810 

1875.1728 

8.65202 

8 

AMS 

135.5817 

395.3242 

1865.3161 

0.00000 

AMS  closed-loop 

146.4793 

392.5993 

1886.3500 

6.55688 
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Tabic  8:  Mode  43  Results 


Case 

Method 

Coast  (sec) 

2nd  Burn  (sec) 

Final  Mass  (kg) 

True  Anomaly  (deg) 

OTIS 

34.1034 

334.7211 

2415.6628 

0.00000 

9 

AMS 

18.1730 

334.3978 

2416.6999 

0.00001 

AMS  closed-loop 

26.7410 

331.8191 

2438.4000 

0.00688 

OTIS 

159.6698 

346.2189 

2311.6721 

0.00000 

10 

AMS 

120.9575 

346.9562 

2303.0468 

0.00000 

AMS  closed-loop 

130.6213 

344.5979 

2320.7500 

0.05877 

OTIS 

307.5751 

390.7985 

1908.4785 

0.00000 

11 

AMS 

219.6417 

390.9169 

1905.2019 

0.00001 

AMS  closed-loop 

238.4151 

390.0967 

1904.4500 

0.33734 

OTIS 

36.9882 

387.9814 

1933.9577 

0.00000 

12 

AMS 

35.0215 

386.7056 

1943.3146 

0.00000 

AMS  closed-loop 

50.0246 

384.6953 

1958.7500 

0.01839 

OTIS 

166.2501 

396.8354 

1853.8786 

0.00000 

13 

AMS 

112.2931 

396.3527 

1856.0081 

0.00000 

AMS  closed-loop 

126.6136 

394.5985 

1868.2500 

0.01075 

OTIS 

293.3618 

431.1413 

1543.6039 

0.00000 

14 

AMS 

187.5651 

430.6763 

1545.3795 

0.00000 

AMS  closed-loop 

206.4970 

430.4641 

1542.4500 

0.06999 

OTIS 

170.3501 

348.4420 

2291.5658 

0.00000 

15 

AMS 

106.7140 

350.5581 

2270.4490 

0.41001 

AMS  closed-loop 

128.8965 

346.3061 

2302.6500 

0.11132 

OTIS 

172.6211 

397.0370 

1852.0556 

0.00000 

16 

AMS 

115.9672 

395.5898 

1862.9119 

0.00000 

AMS  closed-loop 

130.4734 

393.7107 

1877.3000 

0.08649 

OTIS 

27.7985 

302.1463 

2710.2808 

N/A 

17 

AMS 

21.4543 

301.3659 

2715.6387 

N/A 

AMS  closed-loop 

26.5988 

298.4510 

2737.0500 

N/A 

OTIS 

174.1365 

315.1391 

2592.7694 

N/A 

18 

AMS 

132.5598 

315.9119 

2583.9970 

N/A 

AMS  closed-loop 

139.9287 

313.4425 

2601.3000 

N/A 

OTIS 

319.0053 

366.0883 

2131.9660 

N/A 

19 

AMS 

240.3386 

365.9951 

2130.7439 

N/A 

AMS  closed-loop 

258.0000 

365.2223 

2130.7000 

N/A 
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Figure  16:  Altitude  and  velocity  profiles  of  the  AMS  closed-loop  burn-coast-burn  ascent 
trajectories  for  cases  9-11 


Figure  17:  Pitch  and  yaw  angles  along  the  AMS  closed- loop  burn-coast-burn  ascent  trajec¬ 
tories  with  respect  to  the  launch  plumbline  (guidance)  frame  for  cases  9-11 
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600 


Figure  18:  AMS  and  OTIS  altitude  and  velocity  comparison  for  case  2 


Figure  19:  AMS  and  OTIS  altitude  and  velocity  comparison  for  case  3 


72 


^Perigee 


Insertion  Point 


Figure  20:  AMS  closed-loop  ascent  trajectory  and  target  insertion  orbit  for  case  6 


Figure  21:  AMS  closed-loop  ascent  trajectory  and  target  insertion  orbit  illustrating  true 
anomaly  at  insertion  for  case  6 
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11.3  Mode  41  /  51  Comparisons 

Mode  41  consists  of  four  orbital  insertion  conditions  defined  by  the  desired  target  orbital 
element  values  for  the  semi-major  axis  a*,  eccentricity  e*,  inclination  i* ,  and  ascending  node 
Q*.  This  is  a  free  attachment  Mode  in  that  the  final  flight  path  angle  7/  is  not  constrained. 
Mode  51  requires  the  same  orbital  insertion  conditions  as  Mode  41  in  addition  to  constraining 
the  final  flight  path  angle  to  zero  such  that  the  insertion  point  is  the  perigee  of  the  target 
orbit.  Modes  41  and  51  are  similar  to  Modes  31  and  43  respectively  and  differ  only  from 
the  additional  constraint  on  the  ascending  node.  This  fixes  the  target  orbital  plane  leaving 
the  only  unknown  the  location  of  the  perigee  on  the  orbit.  A  set  of  test  cases  with  target 
orbital  insertion  conditions  as  well  as  an  identical  set  for  Mode  51  including  the  final  7/  =  0 
constraint  is  given  below  in  table  9  .  In  addition,  three  circular  target  orbits  for  Mode  51 
are  listed.  This  set  of  test  cases  is  identical  to  those  of  31  and  43  with  the  addition  of  a  fixed 
ascending  node  fl*.  It  should  be  noted  that  for  the  initial  conditions  used  for  all  test  cases 
in  this  work,  target  orbits  with  unconstrained  ascending  nodes  converge  to  approximately 
255°  for  inclinations  of  51.6°  and  185°  for  inclinations  of  28.5°. 


Table  9:  Mode  41  &  51  Target  Orbits 


Case 

Mode 

perigee  alt.  (km) 

e 

i  (deg) 

7/  (deg) 

R  (deg) 

20 

41 

300 

0.1 

51.6 

Free 

250 

21 

41 

500 

0.1 

51.6 

Free 

250 

22 

41 

1000 

0.1 

51.6 

Free 

250 

23 

41 

300 

0.3 

51.6 

Free 

250 

24 

41 

500 

0.3 

51.6 

Free 

250 

25 

41 

1000 

0.3 

51.6 

Free 

250 

26 

41 

500 

0.1 

28.5 

Free 

180 

27 

41 

500 

0.3 

28.5 

Free 

180 

28 

51 

300 

0.1 

51.6 

0.0 

250 

29 

51 

500 

0.1 

51.6 

0.0 

250 

30 

51 

1000 

0.1 

51.6 

0.0 

250 

31 

51 

300 

0.3 

51.6 

0.0 

250 

32 

51 

500 

0.3 

51.6 

0.0 

250 

33 

51 

1000 

0.3 

51.6 

0.0 

250 

34 

51 

500 

0.1 

28.5 

0.0 

180 

35 

51 

500 

0.3 

28.5 

0.0 

180 

36 

51 

1000 

0.0 

51.6 

0.0 

250 

37 

51 

500 

0.0 

51.6 

0.0 

250 

38 

51 

500 

0.0 

51.6 

0.0 

250 

Results  for  Mode  41  and  51  test  cases  in  table  9  are  listed  in  table  10  and  11  respectively. 
Again,  both  the  second  vacuum  burn  time  and  mass  are  listed  for  convenience.  In  terms  of 
vehicle  performance,  all  methods  match  quite  closely  as  was  the  case  for  Modes  31  and  43. 
The  open-loop  AMS  trajectory  continues  to  match  the  full  non-linear  gravity  AMS  closed- 
loop  simulation  which  is  evidence  the  AMS  method  has  no  problem  handling  the  additional 
constraint  on  the  problem.  For  Modes  31  and  43  it  was  seen  that  for  higher  orbital  insertion 
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altitudes,  the  longer  the  coast  arc  tends  to  be,  and  the  more  critical  is  its  length  to  the 
performance  of  the  vehicle.  This  trend  is  again  seen  for  Mode  51  wherein  the  insertion  point 
must  be  the  perigee  of  the  target  orbit.  However,  from  the  results  in  table  10,  longer  coast 
arcs  for  higher  insertion  altitudes  is  not  the  typical  outcome  for  Mode  41.  This  is  explained 
by  the  variation  of  optimal  insertion  points  found  in  terms  of  the  true  anomaly  at  insertion 
in  addition  to  the  unconstrained  target  orbit  perigee  direction  which  can  be  much  different 
from  the  corresponding  Mode  51  cases. 

The  ascending  node  constraint  enforced  in  Modes  41  and  51  appears  to  be  the  dominate 
constraint.  Enforcing  the  ascending  node  constraint  fl  in  addition  to  the  semi-major  axis 
a,  eccentricity  e,  and  inclination  %  fully  defines  the  target  orbital  plane  leaving  the  only  un¬ 
known  the  perigee  direction.  It  can  be  seen  from  table  9  that  the  target  orbit  ascending  node 
for  Modes  41  and  51  are  set  at  approximately  a  5°  offset  from  their  natural  unconstrained 
converged  values  for  the  initial  conditions  used  in  this  work.  Enforcing  this  constraint  has  a 
very  noticeable  impact  on  vehicle  performance.  Every  Mode  41  and  51  case  results  in  a  sig¬ 
nificantly  lower  final  mass  than  the  corresponding  Mode  31  and  43  case  as  shown  in  tables  10 
and  11.  Constraining  the  ascending  node  to  the  values  in  table  9  forces  the  launch  vehicle  to 
make  substantial  out  of  plane  maneuvers  during  ascent  to  align  itself  with  the  target  orbital 
plane.  These  out  of  plane  maneuvers  require  longer  second  burn  times  causing  a  significant 
decrease  in  final  vehicle  mass  at  insertion.  For  Modes  31  and  43  it  was  shown  (see  Fig.  17) 
that  the  yaw  angle  is  generally  very  small  during  the  ascent  due  to  the  initial  launch  azimuth 
direction.  Expensive  out  of  plane  maneuvers  result  in  much  larger  yaw  angles  during  the 
ascent  as  illustrated  in  Figure  23  with  a  comparison  of  case  33  and  corresponding  case  14. 
Figure  22  shows  the  altitude  and  velocity  profiles  for  these  cases.  As  seen  in  the  figures, 
both  cases  have  very  similar  velocity  and  altitude  profiles.  In  case  33,  however,  the  large 
yawing  motion  during  the  ascent  results  in  a  longer  second  burn  time  by  approximately  7 
seconds  translating  into  a  final  mass  reduction  of  approximately  63  kg.  Further  verifying 
these  results,  Figures  24  and  25  show  two  views  of  the  optimized  trajectory  and  target 
insertion  orbit  for  Mode  41  test  case  25.  Again,  easily  seen  from  the  figures,  the  launch  ve¬ 
hicle  performs  substantial  yawing  maneuvers  to  align  and  insert  itself  into  the  desired  target 
orbit.  These  results  are  consistent  for  both  the  AMS  method  and  OTIS. 

Similar  to  Modes  31  and  43,  the  performance  index  on  the  final  mass  in  the  optimal 
burn-coast-burn  problem  continues  to  be  fairly  insensitive  to  the  coast  time.  Significant 
differences  in  the  converged  coast  time  between  OTIS  and  the  AMS  method  does  not  re¬ 
sult  in  significant  differences  in  orbital  insertion  mass.  This  disagreement  in  optimal  coast 
times  is  again  more  noticeable  for  higher  altitude  orbital  insertions.  The  small  changes  in 
insertion  mass  for  disproportionately  large  changes  in  coast  time  is  in  part  due  to  the  “flat 
optimum”  characteristics  of  the  optimal  ascent  problem  as  mentioned  previously.  Further 
understanding  of  the  terminal  Mode  constraints  can  help  explain  this  “flat  optimum”  most 
often  resulting  in  optimal  coast  time  disagreement  between  the  AMS  and  OTIS  methods. 
Is  is  to  be  noted  and  understood  that  the  desired  target  orbits  are  not  fully  defined.  Mode 
41  and  51  does  fix  the  target  orbital  plane  as  mentioned,  but  does  not  constrain  the  perigee 
location.  Given  this,  for  any  test  case  the  AMS  and  OTIS  methods  may  not  insert  into  the 
same  orbit  and  in  general  they  will  not,  however,  both  satisfy  the  terminal  constraints  at 


75 


insertion.  This  freedom  helps  explain  why  both  methods  can  agree  in  terms  of  vehicle  perfor¬ 
mance  even  with  substantial  disagreement  in  optimal  coast  time.  This  can  be  demonstrated 
by  analyzing  a  Mode  51  case  wherein  the  orbital  plane  and  insertion  point  on  the  orbit  are 
fixed.  Figure  26  shows  the  AMS  and  OTIS  trajectories  and  insertion  orbits  for  Mode  51 
case  33.  Both  insertions  must  be  at  the  perigee  of  the  target  orbit,  and  from  table  11,  it  is 
seen  that  the  optimized  coast  arc  for  the  OTIS  trajectory  is  substantially  longer  that  that 
found  by  the  AMS  method.  For  this  reason,  it  is  expected  that  the  OTIS  insertion  point  be 
somewhat  downrange  of  the  AMS  insertion  and  indeed  this  is  the  result.  Figure  27  illustrates 
the  “flat  optimum”  characteristic  of  the  problem,  again  for  case  33,  by  sweeping  the  coast 
time  over  a  200  second  window  bracketing  the  original  converged  coast  times  found  by  both 
the  OTIS  and  AMS  methods.  As  can  be  seen  from  figure,  the  variation  in  final  burn  time 
is  very  insensitive  to  the  coast  duration  in  this  range.  Discussed  above,  the  only  remaining 
unconstrained  orbital  insertion  parameter  is  the  argument  of  perigee,  the  variation  of  which 
is  also  shown  in  the  figure.  Easily  seen  from  the  figure,  increasing  the  coast  time  simply 
pushes  the  perigee  direction  and  insertion  point  further  downrange  from  the  launch  site,  and 
visa  versa. 

In  regards  to  the  solutions  obtained  from  OTIS,  a  few  comments  should  be  noted.  In 
many  cases,  adjusting  the  guessed  final  burn  time  values  was  required  to  obtain  the  optimal 
solution  found.  Further  change  of  these  values  could  prevent  OTIS  from  retrieving  the  same 
solution,  or  in  some  situations  even  a  similar  solution.  It  was  observed  that  adjusting  the 
coast  and  or  final  burn  time  bounds  would  additionally  have  a  significant  impact  on  the 
converged  solution.  This  prevented  obtaining  a  sweep  solution  simply  by  adjusting  the  coast 
time  duration  while  leaving  all  other  parameters  unchanged.  In  addition,  scaling  of  problem 
parameters  highly  influenced  convergence  rates  and  final  solutions.  It  is  recognized  that  OTIS 
is  a  very  general  trajectory  optimization  software  allowing  for  much  user  configuration,  hence 
and  a  more  complete  knowledge  of  its  proper  use  in  regards  to  this  specific  problem  may 
have  relieved  some  of  these  issues. 
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Tabic  10:  Mode  41  Results 


Case 

Method 

Coast  (sec) 

2nd  Burn  (sec) 

Final  Mass  (kg) 

True  Anomaly  (deg) 

OTIS 

259.9649 

351.9201 

2260.1084 

28.88764 

20 

AMS 

219.0324 

352.1562 

2255.9867 

20.38100 

AMS  closed-loop 

231.6762 

349.9753 

2275.5000 

19.04381 

OTIS 

270.9058 

361.7805 

2170.9273 

11.64281 

21 

AMS 

191.2724 

363.3863 

2154.3536 

0.98688 

AMS  closed-loop 

204.5619 

364.9046 

2139.7500 

1.61218 

OTIS 

299.9342 

398.2460 

1841.1204 

8.97000 

22 

AMS 

211.8645 

401.6461 

1808.1027 

10.88184 

AMS  closed-loop 

227.6674 

400.0388 

1813.9500 

12.15001 

OTIS 

265.5381 

399.2454 

1832.0815 

10.57011 

23 

AMS 

210.3717 

397.9093 

1841.9211 

6.07592 

AMS  closed-loop 

227.1023 

396.1889 

1850.1500 

5.03092 

OTIS 

222.2749 

407.0645 

1761.3633 

2.34400 

24 

AMS 

158.6010 

407.3729 

1756.2756 

3.49641 

AMS  closed-loop 

174.9292 

404.9382 

1777.7500 

4.93180 

OTIS 

285.9631 

442.5412 

1440.4989 

11.85391 

25 

AMS 

149.9988 

438.8672 

1471.2519 

11.49603 

AMS  closed-loop 

164.2296 

436.9191 

1488.1500 

14.74207 

OTIS 

130.6578 

345.4328 

2318.7820 

12.48016 

26 

AMS 

88.9099 

346.4386 

2307.7307 

12.04604 

AMS  closed-loop 

94.3723 

343.9121 

2329.8000 

15.26107 

OTIS 

219.2165 

396.1688 

1859.9077 

0.00000 

27 

AMS 

15.7832 

396.0942 

1858.3475 

12.37148 

AMS  closed-loop 

21.0860 

393.6329 

1877.3000 

16.89375 
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Tabic  11:  Mode  51  Results 


Case 

Method 

Coast  (sec) 

2nd  Burn  (sec) 

Final  Mass  (kg) 

True  Anomaly  (deg) 

OTIS 

188.2341 

356.1665 

2221.7026 

0.23231 

28 

AMS 

168.6447 

356.3202 

2218.3019 

0.00000 

AMS  closed-loop 

177.0205 

354.2091 

2230.2500 

0.16114 

OTIS 

228.3949 

366.3698 

2129.4199 

0.22885 

29 

AMS 

188.8833 

363.5015 

2153.3116 

0.00001 

AMS  closed-loop 

198.4591 

361.0573 

2166.9000 

0.01260 

OTIS 

326.8909 

398.5508 

1838.3637 

0.22111 

30 

AMS 

230.6249 

401.1288 

1812.7843 

0.00001 

AMS  closed-loop 

248.3052 

399.9408 

1823.0000 

0.02116 

OTIS 

182.9739 

399.9333 

1825.8598 

0.19379 

31 

AMS 

170.6421 

399.4623 

1827.8661 

0.00000 

AMS  closed-loop 

182.1409 

397.6953 

1841.1000 

0.03067 

OTIS 

235.0993 

405.9083 

1771.8204 

0.14380 

32 

AMS 

180.8982 

405.2135 

1761.5226 

0.00000 

AMS  closed-loop 

193.6338 

404.8655 

1777.7500 

0.00441 

OTIS 

317.0431 

438.8030 

1474.3089 

0.13850 

33 

AMS 

201.6559 

437.9284 

1479.7483 

0.00000 

AMS  closed-loop 

220.8111 

437.2554 

1479.1000 

0.09196 

OTIS 

151.4102 

347.6747 

2298.5052 

0.22882 

34 

AMS 

122.5030 

346.2990 

2308.9943 

0.00000 

AMS  closed-loop 

132.1887 

343.9411 

2329.8000 

0.00545 

OTIS 

158.9689 

397.0465 

1851.9692 

0.14366 

35 

AMS 

114.0376 

395.8322 

1860.7183 

0.00000 

AMS  closed-loop 

128.4066 

394.1371 

1868.2500 

0.00323 

OTIS 

200.4568 

330.3483 

2455.2113 

N/A 

36 

AMS 

171.1778 

330.8833 

2448.5065 

N/A 

AMS  closed-loop 

178.0483 

328.7342 

2465.5500 

N/A 

OTIS 

159.8574 

339.3094 

2374.1639 

N/A 

37 

AMS 

196.6761 

337.2291 

2391.0771 

N/A 

AMS  closed-loop 

204.6825 

334.6980 

2411.2500 

N/A 

OTIS 

337.8690 

380.0354 

2005.8236 

N/A 

38 

AMS 

248.9145 

378.2907 

2019.4688 

N/A 

AMS  closed-loop 

265.7263 

377.1271 

2022.1000 

N/A 
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Figure  22:  Altitude  and  velocity  profiles  of  the  AMS  closed-loop  burn-coast-burn  ascent 
trajectories  for  cases  14  and  33 


Figure  23:  Pitch  and  yaw  angles  along  the  AMS  closed-loop  burn-coast-burn  ascent  trajec¬ 
tories  with  respect  to  the  launch  plumbline  (guidance)  frame  for  cases  14  and  33 
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Figure  24:  View  1:  AMS  closed-loop  ascent  trajectory  and  target  insertion  orbit  for  case  25 
illustrating  large  out  of  plane  motion  (orbit  shading  for  visual  convenience) 
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Figure  25:  View  2:  AMS  closed-loop  ascent  trajectory  and  target  insertion  orbit  for  case  25 
illustrating  large  out  of  plane  motion  (orbit  shading  for  visual  convenience) 
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Figure  26:  AMS  closed-loop  and  OTIS  ascent  trajectories  and  orbital  insertion  perigee  di¬ 
rection  for  case  33 


Coast  Time  (sec) 


Coast  Time  (sec) 

Figure  27:  OTIS  coast  time  sweep  for  case  33 
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11.4  Mode  46  Results 


Mode  46  involves  four  insertion  conditions  defined  by  the  desired  final  position  magnitude 
rjl,  velocity  magnitude  Vf,  inclination  i*7  and  flight  path  angle  'yj.  This  Mode  is  unlike 
the  others  in  that  it  is  not  defined  from  typical  orbital  elements  and  the  conditions  need 
not  represent  values  on  any  target  orbit.  Rather,  this  Mode  provides  the  ability  to  define  a 
desired  launch  abort  energy  condition  where,  if  met,  the  launch  vehicle  can  glide  unpowered 
if  necessary  to  a  determined  landing  site.  Conditions  for  three  cases  are  listed  in  table  12. 


Table  12:  Mode  46  Target  Orbits 


Case 

Mode 

altitude  (km) 

velocity  (m/s) 

i  (deg) 

If  (deg) 

39 

46 

122 

7500 

51.6 

-1 

40 

46 

122 

7500 

51.6 

-2 

41 

46 

100 

7000 

40 

-1 

The  results  for  these  three  cases  are  listed  in  table  13.  All  methods  agree  closely,  both  in 
optimized  coast  length  and  final  burn  time,  however  the  OTIS  trajectories  have  somewhat 
lager  second  burn  times  resulting  in  reduced  final  mass.  The  initial  conditions  used  for  these 
cases  are  the  same  as  all  previous  cases  and  are  conditions  taken  from  a  typical  ascent.  This 
results  in  a  higher  initial  altitude  than  would  most  likely  be  experienced  in  an  abort  scenario, 
thus  increasing  the  convergence  difficulty.  Figure  28  shows  both  the  AMS  and  OTIS  altitude 
and  velocity  profiles  as  well  as  the  flight  path  angle  time  histories  for  case  40. 


Table  13:  Mode  46  Results 


Case 

Method 

Coast  (sec) 

2nd  Burn  (sec) 

Final  Mass  (kg) 

OTIS 

0.9908 

272.5345 

2978.1001 

39 

AMS 

0.0000 

267.9072 

3018.4401 

AMS  closed-loop 

0.1315 

265.6799 

3035.7000 

OTIS 

0.3931 

267.7783 

3021.1169 

40 

AMS 

0.1763 

263.2835 

3060.2844 

AMS  closed-loop 

0.0392 

260.8845 

3080.9500 

OTIS 

0.9037 

226.0943 

3398.1224 

41 

AMS 

0.0007 

222.9545 

3425.2613 

AMS  closed-loop 

0.0318 

221.4846 

3433.9000 
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Figure  28: 
case  40 


AMS  and  OTIS  altitude,  velocity,  and  flight  path  angle  profile  comparison  for 
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12 


Evaluating  the  Complete  Ascent 
Guidance  Algorithm 


In  this  part  the  endo-atmospheric  and  exo-atmospheric  ascent  guidance  algorithms  are  com¬ 
bined  to  generate  complete  optimal  ascent  trajectories  and  guide  the  vehicle  from  liftoff  to 
orbital  insertion.  To  demonstrate  the  applicability  of  the  algorithms  to  both  winged  launch 
vehicles  and  conventional  axisymmetric  launch  vehicles,  a  two-stage  winged  fully  reusable 
launch  vehicle  model  and  the  Ares  I  Crew  Launch  Vehicle  (CLV)  currently  under  develop¬ 
ment  by  NASA  are  used.  The  launch  site  is  at  KSC.  For  testing  purposes  114  different  wind 
profile  pairs  are  used  which  are  based  on  measurement  data  at  KSC  in  February.  These 
wind  profiles  are  generated  by  the  Natural  Environments  Branch  at  NASA  Marshall  Space 
Flight  Center.  Each  profile  provides  altitude-dependent  wind  speed  and  direction.  The 
first  of  each  pair  is  a  smoothed  (filtered)  prohle  of  measured  wind  data;  the  second  is  un¬ 
filtered  and  measured  on  average  2  hours  later.  The  smoothed  prohle  in  each  pair  is  used 
in  the  guidance  solution,  and  the  unfiltered  prohle  is  used  in  simulation  of  the  trajectory. 

In  Monte  Carlo  simulations  the  wind  pair  is  uniformly  randomly  selected  for  a  trajectory  to 
simulate  the  “day-of-launch  wind” .  In  addition,  various  environmental  and  vehicle  modeling 
uncertainties  are  added  in  the  simulations.  But  no  navigation  dispersions  are  considered  in 
the  simulations.  In  each  simulation,  the  guidance  is  given  only  the  nominal  models  of  the 
vehicle  and  atmosphere,  and  the  smoothed  (and  two-hour  earlier)  prohle  of  the  randomly 
selected  “day-of-launch”  wind  prohle.  No  additional  guidance  I-loads  are  adjusted  from  one 
trajectory  to  another,  and  no  pre-launch  guidance  computation  or  update  is  performed. 

12.1  Application  to  Two-Stage  Winged  Reusable  Launch 
Vehicle 

The  vehicle  is  a  two-stage  winged  fully  reusable  launch  vehicle  similar  to  what  is  shown  in 
Fig.  29.  The  hrst  stage  has  4  engines,  and  second  stage  has  one.  Some  of  the  basic  data  for 
the  vehicle  are  given  in  Table  14. 

The  target  orbit  data  are  given  in  Table  15.  No  other  inequality  trajectory  constraints  are 
enforced.  In  addition  to  the  randomly  selected  winds,  other  dispersions  used  in  the  Monte 
Carlo  simulations  are  summarized  in  Table  16. 
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Figure  29:  Two-stage  winged  fully  reusable  launch  vehicle 
Table  14:  Vehicle  data  of  the  two-stage  reusable  launch  vehicle 


Parameter 

First  Stage 

Second  Stage 

Vacuum  thrust 

4  x  2,086,883  N 

990,280.0  N 

Isp 

335.0  sec 

350.0  sec 

Initial  mass 

574,023  kg 

93,725kg 

Empty  mass 

96,807.3  kg 

17,190.8  kg 

Propellant  mass 

382,671.0  kg 

76,534.0  kg 

The  dispersions  listed  in  Table  16  are  some  of  the  representative  dispersions  commonly 
used  in  launch  simulations.  The  main  objective  in  this  section  is  to  evaluate  the  robustness 
of  the  ascent  guidance  algorithm  in  the  inevitable  presence  of  winds  and  these  appreciable 
dispersions.  For  the  algorithm  to  qualify  to  be  considered  for  on-board  applications,  it  must 
be  able  to  ensure  convergence  in  the  face  of  such  dispersions  without  any  knowledge  of  them. 
Also  to  be  demonstrated  is  the  potential  for  automated  closed-loop  ascent  guidance  without 
any  case-dependent  guidance  I-load  update,  which  is  a  key  requirement  for  responsive  launch, 
by  showing  that  the  algorithm  can  correctly  guide  the  vehicle  to  the  desired  orbit  in  every 
simulation  by  simply  “pushing  the  ignition  button” . 

A  total  of  200  Monte  Carlo  simulations  are  shown  in  Figs.  30-36.  Figures  30  and  31 


Table  15:  Target  orbit  for  the  2-stage  RLV  (launch  from  KSC) 
Orbit  Parameter  i  (deg)  e  altitude  (km)  ascending  node  (deg) 


Value 


49.1  0.0 


200.0 


-107.5 


Table  16:  Dispersions  used  in  Monte  Carlo  simulations  for  the  two-stage  RLV  (all  Gaussian 
distributions) 


Parameter 

Dispersion  Mean 

one-sigma  value 

Liftoff  mass 

5,740,234.5  kg 

1,912.2  kg 

Mass  at  first  stage  burnout 

202,379.1  kg 

674.5  kg 

Mass  flow  rate  for  first  stage 

-2,547.6  kg/sec 

8.5  kg/sec 

Vacuum  thrust  for  first  stage 

834,711.6  N 

2,3803.5  N 

Mass  flow  rate  for  second  stage 

-289.2  kg/sec 

0.96  kg/sec 

Thrust  for  second  stage 

990,280.0  N 

2,823.8  N 

Scaling  factor  for  dispersions  in  aero  coefficient  C  a 

0.0 

0.0333333 

Scaling  factor  for  dispersions  in  aero  coefficient  Cn 

0.0 

0.0333333 

depict  the  geodetic  altitude  and  inertial  velocity  along  these  dispersed  trajectories.  From  the 
figures  it  is  evident  that  all  the  trajectories  reach  the  specified  orbit  successfully,  even  though 
we  will  leave  the  examination  of  the  actual  orbital  insertion  statistics  to  the  next  section. 
The  profiles  of  angle  of  attack  profiles  are  plotted  in  Fig.  32.  The  initial  large  variations  of  a 
are  caused  by  the  fact  the  earth-relative  velocity  is  very  small  right  after  liftoff.  Any  winds 
will  cause  large  changes  in  the  direction  of  the  relative  velocity.  The  vehicle  flies  a  heads-up 
configuration  so  after  liftoff  a  is  largely  positive  in  endo- atmospheric  ascent  portion.  The 
pitch  and  yaw  angle  profiles  of  the  vehicle  with  respect  to  the  NDE  frame  are  illustrated  in 
Figs.  33  and  34.  Finally  the  variations  of  the  dynamic  pressure  q  and  the  product  a  —  q  are 
shown  in  Figs.  35  and  36.  Note  that  since  neither  is  constrained,  they  reached  quite  large 
values,  especially  for  aq. 
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Figure  30:  Geodetic  altitude  profiles  along  200  dispersed  RLV  ascent  trajectories 
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Figure  31:  Inertial  velocity  profiles  along  200  dispersed  RLV  ascent  trajectories 
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Figure  32:  Anlgle  of  attack  profiles  along  200  dispersed  RLV  ascent  trajectories 


Figure  33:  Pitch  angle  profiles  along  200  dispersed  RLV  ascent  trajectories 
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Figure  34:  Yaw  angle  profiles  along  200  dispersed  RLV  ascent  trajectories 


Figure  35:  Dynamic  pressure  profiles  along  200  dispersed  RLV  ascent  trajectories 
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Figure  36:  Dynamic  pressure  profiles  along  200  dispersed  RLV  ascent  trajectories 
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12.2  Application  to  Ares  I  Crew  Launch  Vehicle 

12.2.1  The  Ares  I  CLV 


The  Ares  I  CLV  is  a  two-stage  vehicle  currently  under  development  that  will  carry  the  Orion 
Crew  Exploration  Vehicle  (CEV)  to  the  International  Space  Station  (ISS).  Figure  37  shows 
its  conventional,  axisymmetric  conEguration.  The  first  stage  of  the  Ares  I  vehicle  is  a  five- 
segment  solid  rocket  booster  (SRB)  derived  from  the  Shuttle  Solid  Rocket  Booster  (SRB). 
The  vacuum  thrust  profile  and  mass  flow  rate  of  the  SRB  have  a  quite  strong  dependence 
on  time  since  ignition,  and  cannot  be  approximated  by  constants,  as  seen  in  Figs.  38  and  39. 
The  second  stage  engine  is  a  Saturn  J-2  derived  liquid  engine,  dubbed  J-2X.  Its  thrust  and 
mass  flow  rate  are  treated  as  constants  in  guidance.  Other  vehicle  parameters  are  given  in 
Table  17.  The  aerodynamic  coefficients  of  the  vehicle  used  for  guidance  solution  are  smooth 
curve  fittings  to  tabulated  data  as  functions  of  Mach  number  and  angle  of  attack.  The 
vehicle’s  gross  liftoff  weight  is  on  the  order  of  2  million  lbs  with  a  payload  capacity  of  25 
metric  tons. 


Figure  37:  The  Ares  I  Crew  Launch  Vehicle 
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Figure  38:  The  vacuum  thrust  profile  of  the  Solid  Rocket  Booster 


Table  17:  Vehicle  data  for  the  Ares  I  CLV 


Parameter 

First  Stage 

Second  Stage 

Vacuum  thrust 

13,964  kN  at  0.7  sec 

2,088  kN 

Vacuum  Isp 

268.8  sec 

452.1  sec 

Initial  mass 

586,344  kg 

93,725  kg 

Empty  mass 

81,818  kg 

20,422  kg 

Propellant  mass 

504,516  kg 

163,530  kg 

The  International  Space  Station  (ISS)  mission  target  orbit  for  Ares  I  is  a  —11  x  100 
nm  orbit  (with  an  eccentricity  of  0.0159)  at  the  inclination  of  51.6  degrees.  The  insertion 
altitude  into  this  orbit  is  specified  at  70  nm.  The  CEV  will  use  its  own  propulsion  to  reach 
the  ISS  orbit  after  separating  from  Ares  I.  The  orbital  insertion  conditions  are  summarized 
in  the  Table  below 

Table  18:  Ares  I  orbital  insertion  condition  for  ISS  mission  (launch  from  KSC) 


Orbit  Parameter 

i  (deg) 

e  altitude  (nm) 

7  (deg) 

Value 

51.6 

0.0159  70.0 

0.8091 

where  7  is  the  inertial  flight  path  angle  at  orbital  insertion. 

On  the  launch  pad  at  the  Kennedy  Space  Center  (KSC),  the  Ares  I  is  positioned  with 
crew  window  facing  due  east.  It  will  take  about  6  seconds  for  the  vehicle  to  vertically  ascend 
to  clear  the  tower,  so  at  6  seconds  from  liftoff,  the  vehicle  begins  to  maneuver  toward  gravity 
turn.  At  t'2  =  20  seconds  from  liftoff,  the  vehicle  enters  near  gravity  turn.  The  closed-loop 
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Figure  39:  The  vacuum  mass  rate  profile  of  the  Solid  Rocket  Booster 

guidance  update  rate  is  1  HZ  while  the  trajectory  simulation  is  performed  at  100  HZ.  For 
nominal  mission  design  (no  wind)  a  constraint  of 

atq  <  500  (psf-deg)  (1) 

is  imposed,  where  cct  is  the  total  angle  of  attack,  the  angle  between  the  wind-relative  velocity 
and  the  longitudinal  body  axis  of  Ares  I.  In  the  presence  of  winds  and  other  dispersions,  the 
guidance  is  required  to  maintain  the  constraint 

atq  <  3000  (psf-deg)  (2) 

This  constraint  is  critically  important  for  Ares  I  because  of  its  long  slender  configuration 
(about  200  m  in  length).  For  the  structural  integrity  of  Ares  I  during  the  ascent  through 
the  atmosphere,  the  bending  moment  due  to  aerodynamic  forces  must  be  kept  within  design 
bounds.  The  quantity  atq  is  direct  measure  of  the  aerodynamic  bending  moment. 

In  addition  to  the  dispersions  listed  in  Table  19,  altitude-dependent  atmospheric  den¬ 
sity  and  pressure  dispersions  from  the  values  from  the  1976  US  Standard  Atmosphere  are 
generated  from  an  in-house  model  and  used  in  the  simulations  (which  are  unknown  to  the 
guidance).  The  model  can  produce  randomized  atmospheric  dispersions  based  on  location, 
altitude,  and  season.  While  much  simpler  than  NASA’s  GRAM-2007  atmospheric  model 
in  Ref.  40,  this  model  produces  quantitatively  similar  atmospheric  dispersions  as  those  by 
GRAM-2007. 

12.2.2  Nominal  Mission 

All  the  results  presented  here  are  from  closed-loop  3DOF  simulations.  First  the  nominal  as¬ 
cent  without  wind  and  other  dispersions  added  is  examined.  Two  cases  are  compared.  One 
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Table  19:  Dispersions  used  in  Monte  Carlo  simulations  for  the  Ares  I  CLV  (Gaussian  distri¬ 
butions  unless  indicated  otherwise) 


Parameter 

Dispersion  Mean 

one-sigma  value 

Dispersion  scaling  factor  for  mass  flow  rate  for  SRB 

0.0 

0.3333% 

Dispersion  scaling  factor  for  vacuum  thrust  for  SRB 

0.0 

0.3333% 

Dispersion  scaling  factor  for  vacuum  thrust  for  J2-X 

0.0 

0.2667% 

J2-X  specific  impulse 

452.1  sec 

0.73333  sec 

Mass  flow  rate  for  J2-X 

see  above 

see  above 

Scaling  factor  for  dispersions  in  aero  coefficient  Ca 

N/A 

±1% 

Scaling  factor  for  dispersions  in  aero  coefficient  Cn 

N/A 

±1% 

*The  dispersion  scaling  factors  for  Ca  and  are  in  uniform  distributions  within  ±1%. 


is  fully  constrained,  subject  to  (1).  The  other  is  freely  optimized  after  the  vertical  ascent 
at  t\  —  6  sec,  without  any  additional  load  or  path  inequality  constraints.  Figure  40  shows 
the  angle  of  attack  profiles  for  the  two  trajectories.  The  solid  line  is  for  the  constrained 
trajectory.  It  can  seen  that  the  a  profile  is  qualitatively  very  similar  to  that  of  the  a  in  Ref. 
citeDukeman2  which  is  optimized  off-line  and  essentially  the  same  as  the  result  obtained 
from  POST.42  On  the  other  hand,  the  unconstrained  optimal  trajectory  flies  an  a  of  ap¬ 
proximately  5  degrees  (  about  10  times  that  of  the  constrained  one).  The  corresponding 
atq  and  dynamic  pressure  q  during  the  first  stage  are  plotted  in  Fig.  41.  It  is  seen  that  the 
unconstrained  trajectory  experiences  a  peak  atq  of  6000  psf-deg,  12  times  that  along  the 
constrained  trajectory.  The  penalty  for  flying  this  highly  constrained  trajectory  is  about 
0.9%  less  of  mass  delivered,  which  amounts  to  430  kg.  This  is  quite  remarkable,  considering 
the  severity  of  the  constraint.  This  achievement  demonstrates  the  effectiveness  of  the  guid¬ 
ance  optimization  in  the  short  period  of  [£1;  f2]  =  [6,20].  Without  these  initial  maneuvers 
properly  optimized,  one  could  easily  end  up  with  the  constrained  trajectory  suffering  up  to 
15%  —  20%  performance  penalty.  Figure  42  shows  the  comparison  of  the  Euler  angles  of 
the  vehicle’s  body  axes  during  the  ascent.  These  angles  are  with  respect  to  the  North-East- 
Down  (NED)  coordinate  system  at  the  launch  site,  in  the  standard  yaw-pitch-roll  rotation 
sequence.  Not  surprisingly,  the  optimal  yaw  angle  found  is  practically  constant  and  only 
about  2.5  degrees  different  from  the  usual  launch  azimuth 


'i/jaz  =  sin 


-l 


COS  2 
COS  <f>c 


where  i  is  the  target  orbit  inclination  and  (f>c  the  geocentric  latitude  of  the  launch  site.  Figure 
11  illustrates  the  altitude  and  velocity  profiles  of  the  nominal  ascent  with  and  without  the 
atq  constraint  (1). 


12.2.3  Monte  Carlo  Simulations 

Next,  500  dispersed  trajectories  are  simulated  with  the  dispersions  and  randomly  selected 
KSC  winds  described  earlier.  In  each  simulation,  the  guidance  is  given  only  the  nominal 
models  of  the  vehicle  and  atmosphere,  and  the  smoothed  (and  two-hour  earlier)  profile  of 
the  randomly  selected  “day-of-launch”  wind  profile.  No  additional  guidance  I-loads  are 
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10 


Figure  40:  Angle  of  attack  along  nominal  Ares  I  optimal  ascent  trajectories:  top  figure:  the 
complete  a  profile  (stage  1  and  2);  bottom  figure:  a  for  the  first  stage 

adjusted  from  one  trajectory  to  another,  and  no  pre-launch  guidance  computation  or  update 
is  performed.  Upon  ignition,  the  CLV  rieses  for  6  seconds  to  clear  off  the  launch  tower.  At 
this  instant,  the  optimal  closed-loop  guidance  described  in  this  report  takes  over  and  steers 
the  vehicle  through  the  separation  of  the  SRB,  firing  of  the  J2-X  to  the  orbital  insertion. 

The  statistics  on  the  final  (orbital  insertion)  conditions  are  listed  in  Table  20.  Also  in¬ 
cluded  in  the  Table  20  are  the  statistics  of  the  peak  values  of  atq  along  all  these  trajectories. 
The  data  in  the  table  clearly  demonstrate  accurate  orbital  insertion  conditions  and  quite 
consistent  performance  in  terms  of  injected  mass.  Figure  44,  in  which  the  geodetic  altitude, 
inertial  velocity  and  inertial  flight  path  angle  along  the  500  dispersed  trajectories  are  plot¬ 
ted,  confirms  the  successful  orbital  insertion  in  all  cases.  The  high  accuracy  of  the  orbital 
insertion  conditions  seen  in  the  Table  suggests  that  if  navigation  errors  are  included,  the 
orbital  insertion  errors  will  be  on  the  same  order  as  that  of  the  navigation  errors.  This  high 
accuracy,  though,  should  be  viewed  in  the  context  that  no  effects  of  the  uncertainty  caused 
by  “tailoff  thrust”  of  the  J2-X  engine  are  included  when  it  is  commanded  to  shut  off.  Such 
an  uncertainty  can  cause,  for  instance,  the  error  of  a  few  nautical  miles  in  semi-major  axis 
of  the  final  orbit.  Again,  one  of  the  most  critical  concerns  during  atmopheric  ascent  is  the 
bending  moment  in  the  presence  of  the  winds.  The  first  subplot  of  Fig.  45  shows  that  the 
closed-loop  guidance  strategy  is  able  to  keep  the  angle  of  attack  small  throughout  the  flight 
of  the  first  stage.  The  statistics  on  peak  atq  are  reassuring  in  that  the  highest  value  is  about 
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Figure  41:  at-q  and  dynamic  pressure  q  along  the  nominal  Ares  I  optimal  ascent  trajectories 
of  the  first  stage 

2800  psf-deg,  below  the  desired  level  of  3000  psf-deg.  In  fact,  the  data  suggest  that  the 
3-sigma  value  of  the  peak  atq  is  at  about  2300  psf-deg.  A  closer  examination  of  the  data 
reveals  that  the  peak  value  of  2800  psf-deg  occurred  in  only  two  trajectories,  both  with  the 
same  wind  profile.  This  is  a  case  when  the  wind  is  relatively  strong  and  the  measured  wind 
data  used  by  the  guidance  differs  considerably  from  the  wind  encountered  in  ascent.  The 
second  subplot  in  Fig.  45  clearly  shows  the  same  conclusion. 

The  Euler  angles  (pitch,  yaw  and  roll)  of  the  vehicle  body  axes  with  respect  to  the  launch 
NED  frame  during  the  endo-atmospheric  ascent  portion  (which  lasts  about  120  seconds) 
along  the  500  trajectories  are  depicted  in  Fig.  46.  The  complete  pitch  and  yaw  angle  for 
both  stages  of  the  Ares  I  in  all  the  dispersed  trajectories  are  plotted  in  Fig.  47.  Note  that 
the  roll  angle  in  exo-atmospheric  flight  is  immaterial  because  it  does  not  affect  the  direction 
of  the  thrust  vector(thus  roll  angle  is  not  shown  in  Fig.  47),  whereas  in  endo-atmospheric 
flight  the  roll  angle  affects  the  sidelslip  angle,  hence  the  total  angle  of  attack  at. 
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Figure  42:  Euler  angles  along  nominal  Ares  I  optimal  ascent  trajectories 


Table  20:  Statistics  of  500  dispersed  trajectories:  orbital  insertion  conditions  and  peak  atq 
(where  a=  semi-major  axis,  e=eccentricity) 


mean 

std.  deviation 

maximum 

minimum 

r.f  (m) 

6,507,775.52 

0.56 

6,507,776.95 

6,507,773.02 

Vf  (m/s) 

7,797.74 

0.10 

7,798.02271 

7797.34 

If  (deg) 

0.8089 

0.0061 

0.82364 

0.7832 

i  (deg) 

51.60001 

0.00001 

51.60007 

51.59994 

a  (m) 

6,460,810.91 

165.29 

6,461,273.58 

6,460,149.90 

e 

0.01588 

0.0001 

0.01611 

0.01549 

mf  (kg) 

45,827.58 

336.45 

46,588.25 

44,490.092 

peak  atq  (psf-deg) 

1,362.44 

327.91 

2,827.64 

705.82 
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Figure  43:  Altitude  and  velocity  along  nominal  Ares  I  optimal  ascent  trajectories 
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Figure  44:  Geodetic  altitude,  inertial  velocity  and  inertial  flight  path  angle  along  500  dis¬ 
persed  ascent  trajectories  of  Ares  I 
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Figure  45:  Variations  of  first-stage  a  and  atq  along  500  dispersed  ascent  trajectories  of  Ares 
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Figure  46:  Variations  of  first-stage  Euler  angles  along  500  dispersed  ascent  trajectories  of 
Ares  I 
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Figure  47:  Variations  of  pitch  and  yaw  angles  (of  both  stages)  along  500  dispersed  ascent 
trajectories  of  Ares  I 
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Appendix  A 

Costate  Equations  for  3D  Endo- Atmospheric  Optimal  Ascent 

The  notation  used  in  the  following  is  the  same  as  defined  in  Sections  in  Part  ??  of  this 
report.  Let  pr  =  dp/ dr,  Tr  =  dT/dr,  Cp  =  popSrefRo/m(t),  and 

a  _  L/-  Sref  7?( )  Pi )  C A  Pr  ^  ^  S rf  f  1  \\) [){ )  (  '  \  f)r 

pr  =  2Mt)  ’  pr  =  2Mt)  ’ 

CAa  =  dCA/da ,  CNa  =  dCN/da, 
cAMach  =  dCA/dMach,  CNMach  =  dCN/dMach, 

aPvb  =  Pylb  =  PvCOs($  -  a),  apvn  =  Pyln  =  pysin(<f>  -  a) 

Denote  the  altitude-dependent  speed  of  sound  by  Vs  (r).  The  complete  costate  equations  (22) 
and  (23),  after  much  vector  differentiation  and  simplification,  are  given  by 


Pr  =  ~J,PV 


3a 


pvb 


+  apvjj  (  Tr  —  Apr  +  ——CpV^C^ 


*pvn  \  pr 


Nn 


1  r  v2C  — 
2  Vr  P  s  NMach  dr 


2Vr 

r 

r 


■A-Mach 


dr 


“t-  C pidE  X  ]  dpvb 


(< Ca  +  /yV.CA^JV,.  +  \ca.v;-^ 


—  a. 


pvn 


(Cn  +  7rrFvsCNMach)Vr  +  -CNaV; 


214 


2 

da 

&V 


Pv  =  ~Pr  +  Cp 
1 


1  1 

ttpvb(C/ {  +  AMach)  ~  Q"pvn (C/v  +  ^yl ^sCnMoc 


cos  a 

^yCLpvb^Aa  dpvn^ No, )  . 

2  sin  a 


Vr~ 


C'Vr 


p  v  r 


2  sin  a 

where  from  the  first  equation  in  Eq.  (7)  and  Vr  =  V  —  ue  x  r  —  V 


da 


(cosaly. 


V  Vr  sin  a 

Note  that  the  costate  equation  (3)  has  been  simplified  by  recognizing  that 

(hlf  +  lnln  -  h)PV  =  0 


(3) 


{flpvbC ' Aa  QjpvnC No)  lb  (4) 


(5) 


(6) 


because  1&,  ln  and  pv  are  in  the  same  plane  in  the  optimal  solution,  and  1&  and  ln  are  unit 
vectors  that  are  orthogonal  to  each  other. 
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Appendix  B 

Jacobians  of  Thrust  Integrals 


In  the  numerical  solution  of  the  burn-coast-burn  problem  in  this  report,  the  Jacobian  of 
the  constraint  vector  in  Eq.  (34)  can  be  analytically  obtained.  While  much  of  it  just  requires 
standard  (and  careful)  vector  differentiation,  the  most  involved  part  is  perhaps  the  Jacobians 
of  the  thrust  integrals  in  Eqs.  (16)  and  (17).  This  Appendix  provides  the  key  equations  that 
are  very  useful  in  computing  analytically  these  Jacobians.  See  Ref.5  for  some  similar  results. 
Consider  the  thrust  integrals  in  an  interval  [ti,T2]: 

Ic(t 2,ti)  =  f  lPv(()cos(u()AT(()d(  e  R3  (7) 

J  n 

T(t2,ti)  =  f  lPv(()  sin(u() AT(()d(  E  R3  (8) 

J  n 


Let  Ai  =  col(pv(Ti),—pr(ri)/u)  be  the  initial  costate  vector  in  this  interval,  and  h  a  step 
size  parameter  as  follows 

,  T2  -  Ti 

1  ~  N 

where  N  =  2  for  using  the  Simpson’s  rule  to  compute  the  above  integrals  and  N  =  4  for  the 
Milne’s  rule.  Define  the  following  matrices  for  i  —  0, 1, . . . ,  N  in  i?3x6 


A (ih)  =  cos(Luh)/3x3  sin(iu;h)/3x3  (9) 

K(n+ih)  =  tt -  1  ,,  /3x3  -  lpv(Ti  +  ih)ll '(ti  +  ih)  A  (ih)  (10) 

||Pv(ti  +  ih)\\  |_  w  J 

Through  the  thrust  direction  1  Pv(r)  =  Pv(T)/\\Pv(r)\\i  b°fh  Ic  and  Is  are  functions  of 
When  these  integrals  are  evaluated  by  a  numerical  quadrature  method,  it  can  be  shown  that 
the  Jacobians  of  Ic  and  Is  with  respect  to  Ai  are  given  by 

—  '  {  2  —  =  — - —  biAT(ri  +  ih)  cos(ri  +  iuh)K(ri  +  ih)  (11) 

d\i  n,  ^ 

-Air-  —  —  — — — V  biAT(ri  +  ih)  sin(n  +  iuh)K(ji  +  ih)  (12) 

oAi  ns 

s  i=0 

where,  for  the  Milne’s  rule,  we  have 

N  =  4,  ns  =  90,  b0  =  7,  b\  —  32,  b2  =  12,  63  =  32,  b4  =  7 
and  for  the  Simpson’s  rule 

N  =  2,  ns  =  6,  bo  =  1,  b\  =  4,  b2  —  1 
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Appendix  C 

Fixed-Point  Approach  for  Endo- Atmospheric  TPBVP 


The  Fixed-Point  Problem 

In  Section  10,  a  special  form  of  two-point-boundary-value  problem  (TPBVP)  needs  to 
be  solved  for  optimal  endo-atmospheric  ascent.  Define  the  state  vector  x  =  col(r  V )  and 
costate  vector  p  =  col(pr  pv).  Once  the  the  controls  lb  (and  ln)  in  right-hand  sides  of  the 
state  Eqs.  (10),  and  costate  Eqs.  (22)-(23)  are  eliminated  by  the  optimality  conditions  as 
functions  of  x  and  p ,  the  state  and  costate  equations  become  homogeneous.  The  initial  state 
x0  is  given.  Suppose  that  the  burn  out  of  the  first  stage  of  the  launch  vehicle  occurs  at  time 
T\ .  With  a  guessed  optimal  state  at  Ti,  the  exo-atmospheric  algorithm  solves  the  optimal 
vacuum  trajectory  to  the  orbital  insertion  point  and  returns  an  optimal  costate  required  at 
Ti,  denoted  by  p.  Then  the  TPBVP  for  the  endo-atmospheric  portion  must  satisfy  and  the 
special  boundary  condition  p(r\ )  =  p  .  The  special  TPBVP  in  the  endo-atmospheric  portion 
becomes 


x' 

=  f  x(Ti  x,  p) 

(13) 

*(0) 

=  x0 

(14) 

P 

=  fp{r,x,p ) 

(15) 

p(n) 

=  Pin) 

(16) 

where  f  x  and  f  p  are  the  resulting  right-hand  sides  of  the  state  and  costate  equations  once 
the  controls  are  eliminated. 


The  time  interval  [0,  Ti]  is  divided  into  N  sub- intervals  of  the  same  length  h  =  T\/N.  Let 
Xi  =  x{ih )  and  p{  =  p(ih)  be  the  value  of  the  solution  at  the  node  r,  =  ih,  i  —  0, . . . ,  N.  The 
middle  point  between  r*_i  and  r*  is  denoted  by  T;_ i/2  =  Ti  —  h/2.  Therefore  the  differential 
equations  (13)  and  (15)  can  be  approximated  by  central  finite  difference  at  T;_  1/2: 


-(Xi  -  aJi-j)  =  f  x 

\{Pi~Pi-i)  =  fP 
Re-organize  above  equations  as 

Xi  =  Xi-i  +  hf 

Pi- 1  =  Pi  ~  hfP  [ 


P— 1/2: 


Xi  +  Xi- 1  Pi  +  Pi _i 


P—1/2  j 


P-1/2 
P-1/2 


2 

’  2  J 

Xi  +  ®i_! 

Pi+Pi-i\ 

2 

’  2  J 

Xi  +  Xi_ 

1  Pi  +  Pi-i 

2  ’  2 
Xi  +  Xi _!  Pi+Pi_  1 


) 


5 


i  =  1, . . . , N 

i  =  1, . . . ,  N 

i  =  1, . . . , N 

i  =  N,...,  1 


The  boundary  condition  is: 


x0  =  £C(0) 
Pk  =  f>Oi) 


(V) 

(18) 

(19) 

(20) 

(21) 

(22) 
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Let  us  take  a  set  of  initial  guess  on  xpi  =  0,1,...,  N)  and  denote  it  as  xp^(i  = 

0, 1, . . . ,  N ),  where  x(00>  =  x0  as  required  by  the  condition  (21).  Let  pp^(i  =  0, 1, . . . ,  N)  be 

the  corresponding  set  of  costate  vectors  that  satisfy  the  co-state  finite- difference  equations 

Eq.  (20).  By  the  given  boundary  condition  (22)  we  must  have  p ^  =  p(t\).  With  this  known 

condition  and  Xi  replaced  by  xp*\  Eq.  (20)  for  i  =  IV  becomes  a  fixed-point  equation  on 
(o) 

Pn-v 

PN-l=‘nN-l(PN-l)  (23) 

where  »7jv-i(p!v-i)  the  righh  hand  side  of  Eq.  (20)  with  i  =  N  —  1.  Once  the  solution  of 
P/v-i  is  obtained,  p$L2  can  be  found  in  a  similar  way  by  a  fixed-point  equation  in  pffl_2- 
Repeating  this  process  successively,  we  can  conclude  that  after  taking  a  set  of  initial  guess 
on  Xi(i  =  1, . . . ,  N),  the  costate  p\  at  the  i-the  node  can  be  determined  by  a  fixed-point 
equation: 

pT  =  V i(p[0)),  i  =  N-  1, . . . ,  0  (24) 

The  above  process  is  used  to  find  the  solution  for  each  pf  \  %  —  1, . . . ,  N  —  1.  Now  denote 
the  solution  of  the  state  finite-difference  equation  (19)  by  xp-\i  =  0, 1, ...  ,1V)  with  p{  in 
Eq.  (19)  replaced  by  the  just  found  pd°),  where  and  note  that  x(()l>  =  x0  again  as  required. 
In  a  similar  fashion  Eq.  (19)  for  i  —  1  becomes  a  fixed-point  equation  on  cc^: 

*i1}  =  Ci(^i1})  (25) 

where  (^(cc^)  is  the  right  hand  side  of  Eq.  (19)  at  i  =  1. 

After  the  solution  of  x^  is  obtained,  the  fixed-point  equation  for  x:p  is  formed  in  the 
same  way.  Repeating  this  process  successively,  the  state  x(-l>  at  ith  node  can  be  obtained  by 
a  fixed-point  equation: 

*i1)  =  C<(*?)),  i  =  l  ,...,1V  (26) 

Now  we  compare  with  the  initial  guess  x.p\i  =  1, .. ..  .V).  If  they  are  sufficiently 
close  to  each  other,  stop  and  an  endo-atmospheric  ascent  trajectory  satisfying  the  finite- 
difference  equation  system  Eqs.  (19-22)  has  been  found.  Otherwise  using  the  newly  obtained 
xp^ii  =  1, . . . ,  N)  to  replace  xp\i  =  1, . . . ,  N)  in  Eq.  (20),  we  will  repeat  the  procedure  to 
find  pp^ii  =  0, 1, . . . ,  N ).  Therefore  there  is  a  total  of  three  fixed-point  iterations  used  for 
solving  the  formulated  TPBVP:  two  inner-loop  fixed-point  iterations  described  in  Eqs.  (24) 
and  (26)  above,  and  one  outer-loop  fixed-point  iteration  that  determines  the  junction  point 
T\  between  the  endo-  and  exo- atmospheric  portion  of  the  complete  ascent  trajectory. 

Described  above  is  a  brief  review  of  the  fixed-point  formulation  for  solving  the  special 
form  of  TPBVP  in  endo-atmospheric  ascent  portion.  The  convergence  of  the  fixed-point 
iterations  is  still  an  open  question.  In  the  following  some  additional  ideas  are  discussed  to 
facilitate  the  convergence  of  the  fixed-point  iterations. 

Contraction  Mapping 

There  are  totally  three  fixed-point  iterations  in  the  presented  algorithm.  To  avoid  possible 
confusion  and  for  convenience  of  expression,  a  general  fixed-point  iteration  expression  will 
be  used  in  discussion  in  this  section. 
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Consider  a  fixed-point  equation 


*  =  /(*) 

(27) 

The  fixed-point  iteration  applied  to  this  equation  is: 

')  -  f(zu)).  J  =  0,1,2,... 

(28) 

where 

zu)  =  zb\  j  =  1,2,... 

(29) 

The  convergence  of  this  iteration  depends  on  whether  the  mapping  on  the  right  hand  of 
Eq.  (27)  is  a  contractive.  For  a  fixed-point  iteration,  it  is  possible  to  facilitate  the  convergence 
of  the  iteration  through  redefining  the  fixed-point  iteration  with  an  appropriate  coefficient  9 
in  the  following  modified  equation  (cf.  Isaacson  et  al  1994). 

zU+i)  =  0£(j+i)  +  (l  -  9)z^\  j  =  0, 1,  2, . . .  (30) 

With  Eq.  (28)  and  Eq.  (30),  we  have  just  defined  a  new  fixed-point  equation: 

z  =  F(z)  :=9f(z)  +  (l-9)z  (31) 

Evidently  the  solutions  to  Eq.  (31)  and  Eq.  (27)  are  the  same  for  any  9^0.  The 
difference  is  that  in  the  cases  where  fixed-point  iteration  (27)  does  not  converge,  we  may 
be  able  to  find  a  proper  9  to  enable  the  convergence  of  the  fixed-point  iteration  (34).  In 
other  words,  we  are  trying  to  formulate  a  contraction  mapping  z  =  F(z)  instead  of  solving 
the  fixed-point  iteration  (27)  directly.  It  should  be  mentioned  that  9  can  be  either  a  scalar 
or  a  matrix.  See  Ref.37  for  more.  Reference38  proposes  and  compares  several  numerical 
techniques  toward  the  solution  of  the  above  special  TPBVP. 
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